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The  transmission  of  digital  imagery  has  become  a 
necessity  in  recent  times.  Systems  such  as  communications 
and  weather  satellites,  facsimile,  remote  control,  and 
machine  intelligence  can  and  do  make  use  of  data 
compression  techniques  to  reduce  bandwidth  and  power 
consumption. 

Research  on  these  techniques  has  led  to  one  form  of 
image  data  compression  which  achieves  good  image  quality 
for  intraframe  coding  at  low  data  rates.  This  technique 
is  known  as  Adaptive  Stochastic  Picture  Coding  (ASPC) 
which  consists  of  a  one-dimensional  unitary  transform  for 
column-wise  decorrelation  used  in  conjunction  with 
Adaptive  Dif feriential  Predictive  Coding  Modulation  (ADPCM) 
for  the  row-wise  decorrelation,  followed  by  quantization 
to  give  the  desired  data  compression.  This  system 
requires  use  of  quantization  techniques  which  limit  system 
performance.  Optimization  of  adaptive  scalar  quantizers 
and  use  of  vector  quantizers  aid  in  the  adaptation  of  the 
system  to  variations  in  the  image  statistics.  This  report 
represents  a  study  of  such  quantizers  in  the  ASPC  system. 
By  examining  these  quantization  methods,  it  will  be  shown 
that  it  is  vital  that  the  proper  quantizer  be  incorporated 
into  the  system  to  achieve  a  particular  data  rate  at 
desired  distortion  levels. 

r 


LlSLt  Ql  Figures 


1.  Classes  of  Quantizers  .  5 

2.  Uniform  Quantizer  .  7 

3.  Scaled  Uniform  Quantizer  .  8 

4.  Quantizer  Types  . 10 

5.  Companding . 13 

6.  DPCM . 15 

7.  Adaptive  Quantizer  .  17 

8.  General  Quantizer  . 18 

9.  Block  Quantizer  . 24 

10.  Transform  Coding  System  .  25 

11.  Lattice  Quantizer  .  30 

12.  Intra-frame  DPCM . 32 

13.  DPCM . 36 

14.  2~Dim  Transform  Coding . 3  9 

15.  Hybrid  Transf orm/DPCM  Coding . 41 

16.  Optimum  Quantization  .  51 

17.  2-Dim  Vector  Quantizer  .  54 

18.  AS PC  Transmitter  with  Forward  Adaptive  Max  ....  63 

19.  AS  PC  Receiver  with  Forward  Adaptive  Max  . 64 

20.  AS PC  Transmitter  with  Backward  Adaptive  Max  ...  72 

21.  AS  PC  Receiver  with  Backward  Adaptive  Max . 73 

22.  Vector  ASPC  Transmitter  . 75 

23.  Vector  ASPC  Receiver . 76 

24.  ASPC  SNR  Table  .  .  . . 87 

25.  ASPC  Predictor  Performance . 88 

26.  ASPC  Quantizer  Performance . 91 


CHAPTER  ONE 


INTRODUCTION 

The  purpose  of  the  work  reported  here  is  to  examine 
various  quantization  and  coding  methods  for  suitability  as 
an  integral  element  in  an  image  data  compression  and 
transmission  system  —  specifically  the  Adaptive 
Stochastic  Picture  Coding  (ASPC)  system  (sometimes  also 
known  as  the  Adaptive  Hybrid  Picture  Coding  system  - 
AHPC) .  Pursuant  to  this  aim,  it  was  necessary  to  perform 
a  qualitative  analysis  of  system  requirements  and 
performance  objectives.  On  the  completion  of  this 
analysis,  the  results  constituted  a  set  of  guidelines  for 
the  desired  system  operation. 

The  main  criteria  specified  is  the  ability  to 
transmit  a  picture  with  a  good  overall  subjective 
reconstruction  at  low  data  rates.  This  is  rephrased  as  a 
desire  to  transmit  pictures  at  low  data  rates  with 
an  "acceptable"  amount  of  distortion.  As  is  expected, 
built  in  limitations  of  the  system  (ASPC  in  this  case) 
prevent  the  achievement  of  the  theoretical  optimal 
performance.  One  can  however,  "fix"  most  of  the  system 
configuration  to  a  set  standard  and  study  the  effects  of 


replacing  an  element  of  the  system  (here  the  quantizer) 
with  various  experimental  schemes.  The  observed  results 
of  the  system's  performance  for  each  scheme  are  then 
evaluated  in  terms  of  meeting  objective  criteria. 

Since  the  target  of  this  report  was  the  study  of 
the  effects  of  various  quantizers  on  the  AHPC  system,  the 
techniques  of  fundamental  quantizer  design  were  applied, 
and  the  different  configurations  of  quantizers  were 
incorporated  into  the  system.  Extensive  computer 
simulation  of  each  quantizer  implementation  was  made  in 
order  that  the  system  performance  could  be  observed  for  a 
cross-comparison  of  quantizer  schemes. 

Since  very  little  of  the  visual  process  is 
understood  to  date,  many  image  processing  systems  are 
modeled  after  communications  and  speech  processing 
systems.  Although  image  and  speech  systems  exemplify  two 
different  signal  processing  applications,  they  both 
exhibit  similarities  in  the  type  of  signals  they  must  deal 
with  and  the  ways  in  which  these  signals  are  dealt.  Due 
to  the  similarities  of  the  design  of  ASPC  to  speech 
processing  systems,  the  analysis  and  design  of  several  of 
the  quantization  methods  were  paralleled  to  those  used  in 
speech  processing  applications. 

The  following  chapters  will  describe  the  theory  and 
outcomes  of  this  quantizer  study,  beginning  with  the  basic 
fundamentals  of  quantization  theory  appearing  in  Chapter 
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Two,  followed  by  a  description  of  the  ASPC  system  in 
Chapter  Three.  Chapter  Four  presents  the  derivation  of 
the  Max-Lloyd  concept  for  scalar  quantization.  Chapter 
Five  gives  the  derivation  and  design  algorithm  for  vector 
quantization.  Chapter  Six  gives  detailed  explanations  of 
the  various  system  configurations  and  program  simulation 
descriptions.  Chapter  Seven  summarizes  the  results  and 
conclusions  of  the  report  including  performance 
evaluations,  data  rate  calculations,  and  thoughts  about 
future  research  projects  incorporating  these  findings. 
Appendix  A  gives  the  supporting  photographs  of  the  image 
reconstruction  as  well  as  all  error  histograms  between  the 
original  image  and  received  image  as  support  for 
conclusions  drawn  in  Chapter  Seven.  Appendix  B  gives  the 
program  flow  charts  and  listings. 


CHAPTER  TWO 


OVERVIEW  0£ 


THEORY 


In  this  chapter,  the  problems  ~  selecting  the 
appropriate  methods  of  signal  digit  ition  needed  to 
approximate  a  source  signal  waveform  s  addressed.  In 
particular,  an  examination  will  be  made  •'  various  types  of 
signal  quantization  including  theory,  implementation, 
performance,  and  general  applications  of  quantizers  with  the 
goal  of  guiding  the  selection  of  a  quantization/coding 
scheme  for  image  transmission. 

Signal  quantization  may  be  defined  as  a  mapping  of 
samples  from  an  analog  source  signal  into  a  finite  set  of 
values  (constrained  alphabet)  representing  the  samples  in  an 
attempt  to  approximate  or  identify  the  source.  Methods  of 
obtaining  this  mapping  are  of  primary  concern.  It  will  be 
shown  that  the  specific  application,  as  well  as  performance 
and  implementation  criteria,  will  be  of  fundamental 
importance  to  the  selection  of  appropriate  methods. 

Signal  or  waveform  quantizers  have  been  proposed  in 
multitudinous  variety.  However  each  exhibit  certain 
characteristics  which  can  be  classified.  Some  basic 


classifications  of  many  quantizers  are  shown  in  figure  1. 
Each  class  of  quantizers  contains  its  own  particular 
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advantages  and  disadvantages  due  to  theory  and  design 
tradeoffs. 


I.  TIME- INVARIANT  QUANTIZERS 


A.  UNIFORM  QUANTIZATION 


Gersho  [11]  describes  memoryless  N-point  symmetrical 

quantizers  by  specifying  a  set  of  N+l  decision  levels  {x  }, 

k 

k=0,...,N,  and  a  set  of  N  output  points  {y  },  j=l,...N. 

j  til 

When  the  value,  x,  of  an  input  sample  lies  in  the  i 

quantizing  interval,  namely  R  ={x  <x<x  },  the  quantizer 

i  i-1  i 

produces  the  output  value  y  .  Since  y  is  used  to 

i  i 

approximate  samples  in  R  ,  it  follows  that  y <  R  .  The  outer 

i  i  i 

levels  are  chosen  so  that  x<(x  ,x  ).  See  figure  2.  The 

0  N 

description  of  this  class  of  quantizers  certainly  fits  the 
general  definition  stated  previously. 

Probably  the  first  widely-used  quantizers  were  those  of 
the  type  described  by  Gersho,  particularly  the  one¬ 
dimensional  time-invariant  class.  These  provide  easy 

implementation  due  to  their  simplistic  design.  A  quantizer 
of  this  sort  merely  samples  the  signal  to  be  quantized 

(normally  at  a  fixed  uniform  sample  rate)  and  immediately 
produces  a  representative  value  for  each  sample  from  a  fixed 
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set  of  finite,  uniformly  spaced  values  symmetrically  placed 
about  the  average  input  value.  A  typical  input/output  graph 
is  illustrated  in  figure  3. 

Two  basic  styles  of  design  now  surface:  the  mid-tread 
and  the  mid-riser.  The  mid-tread  produces  a  zero  output  for 
signals  at  the  mean  input  level.  The  mid-riser  produces 
output  values  shifting  from  a  negative  value  to  a  positive 
value  at  the  mean  input  level.  See  figure  4.  Note  that 
the  terms  positive  and  negative  could  be  with  respect  to  a 
mean  value  instead  of  zero. 

As  is  intuitively  obvious,  this  type  of  quantizer  would 
perform  best  for  stationary  uniformly  distributed  source 
signals.  If  the  source  statistics  were  to  change  with  time, 
the  quantizer  could  overload  —  producing  useless  output. 
If  we  express  the  quantization  distortion  as  d(x)  =  Q(x)-x 
we  find  the  minimum  distortion  for  this  system  is  produced 
when  the  input  signal,  x,  is  uniformly  distributed.  This  is 
true  regardless  of  whether  it  is  of  the  mid-tread  or  mid¬ 
riser  type. 

B.  NON-UNIFORM  QUANTIZATION 

Most  signals  of  interest  are  seldom  of  the  uniformly 
distributed  class.  Therefore,  to  accommodate  arbitrary 
distributions,  a  design  procedure  must  be  developed. 
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Assuming  the  signals  to  be  stationary,  the  case  of  signals 
with  known  probability  distributions  is  investigated. 
Applications  where  a  specific  probability  density  function 
is  known  to  describe  the  source  signal  can  use  quantizers 
which  are  designed  to  take  advantage  of  this  special 
knowledge.  Derivations  of  quantizers  for  these  signals  are 
presented  in  the  classic  works  by  Lloyd  [24]  and  Max  [251 
and  will  be  discussed  extensively  in  a  later  chapter.  These 
papers  dealt  with  producing  optimum  memoryless  signals  with 
probability  distributions  known  a  priori. 

Quantizers  which  minimized  distortion  defined  by  the 

mean  squared  quantization  error  were  specifically  derived, 

yet  the  same  general  procedure  could  be  used  for  distortion 

measures  of  other  types.  Roe  [33] ,  for  example,  derived 

e 

quantizers  for  all  distortions  of  the  form:  E{error  }. 
Tabulated  results  for  the  decision  levels  and  output  points 
were  produced  for  a  standard  normal  distribution  from  the 
iterative  procedure  derived  [24]  [25].  Others  have  since 

formed  tables  for  other  common  distributions,  such  as  the 
Laplacian  and  Rayleigh,  useful  in  speech,  optical 
holography,  and  image  processing  applications  [22]  [23], 


C.  COMPANDING 

Another  method  for  digitizing  signals  of  this  nature 

has  been  suggested  by  Bennett  [2] .  The  non-uniform 

quantizer  is  modeled  as  a  zero-memory  nonlinearity  function, 

F(x) ,  followed  by  a  uniform  quantizer,  in  turn  followed  by 

-1 

the  inverse  of  the  function,  F  (x) .  The  procedure  is 

called  "companding".  A  diagram  appears  in  figure  5.  The 

first  nonlinear  function,  F(x),  compresses  the  input  signal 
by  spreading  out  the  low  amplitude  portions  and  shrinking 
those  of  high  amplitude.  Since  the  low  amplitudes  have  a 
higher  probability,  this  allows  a  larger  proportion  of 

decision  levels  of  the  uniform  quantizer  to  be  devoted  to 
higher  probability  of  occurence  regions,  and  fewer  to  those 
of  lower  probability  —  increasing  the  accuracy  of  the 
representation.  The  inversion  at  the  end  reverses  the 
process,  producing  the  estimate  of  the  signal.  This 

compression  and  expansion  of  the  signal,  hence  the  name 
"companding",  combine  to  yield  results  comparable  to  those 
of  the  non-uniform  quantizer. 

Holzwarth  and  Smith  [11]  have  made  a  study  of  speech 
transmission  systems  with  widely  varying  power  levels  using 
companding.  They  determined  a  compressor  function  of: 


log  (1- fix  (v) ) 
F  ( x)  =  V - 

log (1+m) 


Since  known  as  a  nfl- law  curve". 


Turning  now  to  the  non-stationary  sources,  it  is  only 
natural  that  research  continue  from  achievements  made  thus 
far.  Certain  applications  such  as  speech  digitization  and 
image  coding  systems  must  deal  with  this  type  of  source. 
Quantizers  which  can  adapt  to  time-varying  source  statistics 
and  still  yield  satisfactory  results  are  vital  to  such  data 
compression  schemes.  Whichever  scheme  one  uses,  such  as 
APC,  DPCM,  or  PPDPCM,  the  quantizer  will  determine  the 
overall  system  performance.  One  common  performance  measure 
is  the  signal-to-quantization  noise  ratio. 

Most  of  the  recent  work  in  speech  digitization  and 
image  coding  systems  concentrate  on  the  use  of  designs 
related  to  the  DPCM  (Dif feriential-Pulse-Coded-Modulation) 
configuration.  Gibson  [15]  gives  a  good  tutorial  on  such 
designs.  The  aspects  of  fixed  and  adaptive  predictors 
utilized  in  such  systems  will  not  be  presented  here,  but  it 
should  be  noted  that  they  too  will  play  an  important  role  in 
system  performance  and  efficiency  by  reducing  signal 
redundancy  such  as  speaker  pitch  in  speech  systems.  The 
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C A  J  DPCM  TRANSMITTER  LOOP 


general  DPCM  differential  encoder  and  decoder  block  diagrams 
are  illustrated  in  figure  6.  Systems  using  these  schemes 
without  any  adaptive  consideration  are  optimal  only  for 
stationary  processes. 

For  the  quantizer  to  adapt  to  the  signal's  statistical 
variations,  it  must  acquire  some  form  of  estimate  of  these 
statistical  variations,  usually  the  signal  variance,  upon 
which  it  is  to  base  its  adaptation.  Two  common 
classifications  of  adaptive  quantizers  are  provided  by 
separating  quantizers  as  to  whether  they  are  forward 
adaptive  or  backward  adaptive.  Forward  adaptive  quantizers 
extract  the  information  for  their  estimate  from  the 
quantizer  input  and  must  transmit  some  additional 
information  regarding  the  subsequent  step  size  variation  to 
the  receiver.  Backward  adaptive  quantizers  utilize  the 
quantized  signal  alone,  eliminating  the  need  for  additional 
transmission.  The  backward  adaptive  quantizer  must  however, 
base  its  estimate  on  cruder  information  due  the  quantization 
noise.  See  figure  7. 

The  basic  idea  is  to  modify  the  quantizer  input  and 

output  step  sizes  in  an  effort  to  match  the  changing  signal 

using  a  signal  variance  estimate.  Goodman  and  Gersho  [17] 

describe  well  the  adaptation.  Again  the  quantizer  of  2N 

output  levels  is  characterized  as  in  figure  8.  It  contains 

a  set  of  N+l  decision  levels  (A(k)^  },  0<i<N,  and  a  set  of  N 

i 
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(AJ  FORWARD  ADAPTIVE  QUANTIZER 


IB)  BACKWARD  ADAPTIVE  QUANTIZER 


Adaptive  Quantizer:  a)  Forward  Adaptive 

b)  Backward  Adaptive 

Figure  7 


General  Quantizer 
Figure  8 


quantization  levels,  {A(k)rj  },  l<i<N.  For  an  input 

i 

z(k),  the  quantizer  is  given  by: 


Q  ( z  (k) )  =  A(k  )rj  z>0;  0<A(k)£  <z(k)<A(k)£  (la) 

i  i-1  i 

and  by  symmetry : 


Q(z(k) )  =  -Q(-z (k) )  z<0.  (lb) 


The  parameters  }  and  {  >?  }  are  forced  to  approximate 

i  i 

the  shape  of  the  desired  probability  density  determined 
through  the  distortion  minimization  analysis  such  as  that  by 
Lloyd  [24]  and  Max  [25]  for  stationary  sources.  The  scaling 
factor,  or  step  size  (k)>0,  which  determines  the  dynamic 
range,  will  be  made  to  vary  with  time.  Several  good 
algorithms  for  determining  step  size  have  been  formulated, 
such  as  these  three  given  by  Stroh  [15] : 


AdO 

A(k) 

AOO 


=  <\ll/k  £  e*  ( k-i) ' , 

•  i-1  2 

=  a  [  (1-a  )  £a  e 
1  2  ‘‘1 2 
1-  a,  £  i"1 

az  1 


1/2 

(k-i) ] 
e(k-i)  , 


(2) 

(3) 

(4) 


for  both  forward  and  backward  adaptive  quantizers  to  be  used 
in  a  DPCM  scheme.  Noll  [15]  proposed  two: 
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J  < 


A(k)  =  a max{ | s (k-i) -s (k-i-1 )  }, 


(5) 


for  first  order  predictors  and 

£  «  21/2 

A(k)  =  a  {  £[s(i)-  £  a  s(i-j)]  }  (6) 

i--l  j*i  j 

where  {a  ,  j=l,...,N}  are  N  predictor  coefficients  for  an 

th  j 

N  -  order  predictor.  Jayant  and  Cummiskey  [15]  obtained 
the  popular  "one-word  memory"  quantizer  with  step  size: 

A(k+1)  =  M(|  Q(k)  | )  A  (k)  (7) 

which  is  a  backward  adaptive  quantizer.  M(>)  is  a  time- 
invariant  multiplier  function  dependent  upon  the  quantized 
signal  as  defined  in  (1)  with  z(k)  =  e(k) ;  e(k)  some 
residual  error  signal  such  as  in  DPCM. 


It  has  been  found  that  an  adaptive  quantizer  based  on  a 
step  size  as  in  (7)  performs  well  for  ideal  channels,  but  is 
quite  susceptible  to  transmission  errors.  Gibson  and 
Wilkenson  [15]  have  shown  that  the  difficulty  can  be 
mitigated  with  the  algorithm  modified  to: 


A(k+1) 


.  .  0 

M<  |Q(k)  |)  •  A  (k) 


(8) 


where  0<  /3<1  will  provide  a  degree  of  robustness.  Gibson 
[15]  discusses  a  closely  related  backward  adaptive  quantizer 
proposed  by  Cohn  and  Melsa  [9]  [10]  where  the  step  size 
adjusts  in  the  same  fashion  as  the  adaptive  Jayant  quantizer 


except  for  two  rarely  used  outer  levels.  The  purpose  of 
which  is  to  allow  rapid  expansion  when  a  pitch  pulse  occurs 
in  a  speech  encoding  system.  It  is  thus  termed  a  pitch- 
compensating  quantizer  (PCQ) ,  and  has  seen  principal 
applications  with  fixed  and  backward  adaptive  predictors. 


B.  CODING 


Once  the  signal  is  quantized,  it  is  usually  encoded  in 
some  fashion  for  transmission.  The  most  common  form  of 
coding  for  an  N-level  quantizer  is  binary  coding.  The  idea 
being  to  assign  a  code  word  of  length  log2N  for  each  output 
level.  This  is  based  on  the  assumptions  that  the  symbols 
are  independent  and  all  quantization  levels  are  equally 
likely.  The  assumption  of  source  symbol  independence  is 
generally  untrue  for  most  applications.  Huang  and 
Schultheiss  [20]  offer  a  technique  to  remedy  this  problem 
and  will  be  discussed  momentarily. 

The  assumption  of  equally  likely  quantization  levels  is 
also  untrue  for  applications  such  as  speech.  Researchers 
have  taken  advantage  of  this  fact  by  designing  variable- 
length  codes  for  the  output  levels.  Using  a  priori 
knowledge,  one  can  assign  code  words  of  short  length  for 
highly  probable  quantization  levels  and  code  words  of  longer 
length  for  those  levels  with  lower  probability.  Variable 


length  coding  is  also  known  as  "entropy  coding",  since  the 
average  code  word  length  can  approach  the  entropy  of  the 
symbols  transmitted  when  the  symbols  are  indeed  independent. 
Hence,  the  entropy  given  as: 


(9) 


becomes  a  lower  bound  for  the  average  code  word  length. 
Berger  [4]  defines  minimum  entropy  quantizers  in  a  rate 
distortion  sense  by  defining  a  rate: 

n 

R  =  -  log  p  (10) 

&  i  2  i 

and  distortion: 


j  H  r  ^  i  j 

d  -  e|y-x|  /  |x-y  |  , 


dF  (x) 

where  {a  }  is  the  set  of  quantization  thresholds  and 


(11) 
{ y  } 


i  i 

the  set  of  reconstruction  levels.  The  optimum  entropy 

quantizer  is  therefore  one  which  minimizes  the  distortion  D 

for  a  fixed  rate  R. 

Berger  [4]  also  discusses  permutation  coding  as 
developed  by  Slepian  [34]  and  introduces  a  theorem  showing 
variable-length  coding  of  quantizer  output  and  permutation 
coding  are  equivalent  in  the  sense  that  their  optimum  R 
versus  D  performances  are  identical. 
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C.  VECTOR 


I 

y 
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All  that  has  been  discussed  previously  is  now  extended 
to  the  next  level  of  sophistication.  This  new  level  is  one 
which  will  deal  with  quantizing  multiple  or  "blocks"  of 
samples  called  quantization  with  memory.  Block  quantization 
is  exactly  that,  quantizing  the  samples  in  blocks.  One  can 
think  of  the  quantizers  discussed  previously  as  one¬ 
dimensional  quantizers,  and  block  quantizers  operating  on 
blocks  of  K  samples  as  K-dimensional .  The  procedure  uses 
blocks  of  quantizers  for  transmission  as  in  figure  9. 

Huang  and  .Schultheiss  (20]  incorporate  this  method  in  a 
system  for  the  transmission  of  correlated  gaussian  random 
variables.  They  use  a  linear  transformation  (P)  to  first 
convert  K  dependent  random  variables  into  K  independent 
random  variables.  These  are  quantized  one  at  a  time  and  the 
output  of  each  quantizer  is  transmitted  by  a  binary  code. 
Lastly,  a  second  KxK  linear  transformation  (R)  constructs 
from  the  quantized  values  the  best  estimate  (in  the  mean- 
squared-error  sense)  of  the  original  values.  They  also  show 
that  the  best  transformation  R  would  be  the  inverse  of  the 
first  transformation  P  The  best  transformation  P  is  the 


transpose  of  the  orthogonal  matrix  which  diagonalizes  the 
moment  matrix  of  the  original  (correlated)  random  variables. 
A  diagram  of  this  process  is  illustrated  in  figure  10. 


SET  OF  K 


Although  the  terms  block  and  vector  quantization  are 
often  used  interchangeably,  when  one  sees  the  term  vector 
quantization  it  is  usually  in  reference  to  a  method  more 
subtle  than  the  brute  force  method  of  K  quantizers.  Gersho 
[13]  states  that  block  or  vector  quantization  is 
intrinsically  superior  to  predictive  coding,  transform 
coding,  and  other  suboptimal  procedures  since  it  achieves 
optimal  rate  distortion  performance  subject  only  to  memory 
(block)  length  of  the  observable  segment  being  encoded.  One 
can  think  of  vector  quantization  as  a  method  of  representing 
a  sequence  of  K  samples  of  the  signal  as  a  vector  in  K- 
dimensional  space.  This  space  is  partitioned  into  a  finite 
number  of  regions,  the  geometry  of  which  is  determined  in  a 
manner  ~ such  that  certain  optimality  criteria  are  achieved. 
Each  region  of  space  is  then  assigned  a  specific  code  word. 
As  the  region  of  space  for  each  sample  vector  is  identified, 
its  representative  code  word  is  selected  and  after  trans¬ 
mission  can  be  decoded  via  a  code  lookup  table  or  "code 
book"  into  an  output  vector  representing  the  original  input 
vector.  Thus  vector  quantization  is  a  surjective  mapping  of 
K-dimensional  input  vectors  to  corresponding  representative 
vectors, 

k  k 

Q  :  R  - ►  Y,  YcR  (12) 
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D.  RANDOM 


AND  CODE  ROOK  DESIGN 


Gersho  [13]  analyzes  the  approach  to  vector 

quantization  by  Linde,  Buzo,  and  Gray  [23]  with  random 

quantizers.  Since  the  only  effective  method  for  the  design 

of  multidimensional  quantizers  is  the  use  of  a  clustering 

algorithm,  their  design  algorithm  uses  a  training  set  of 

random  vectors  generated  from  the  source.  The  training  set 

is  the  collection  [x  ,...,x  }  of  M  independent  observations 

1  M 

of  the  continuously  distributed  random  variable  X  with  a 

specified  joint  density  function  where  M  >>  N  for  an  N  level 

or  N  regional  quantizer.  The  set  of  output  vectors 

[y  ....,y  }  is  determined  by: 

1  H 


y  x  s  (x  ) 

•  i  i  i 

i*l.  J 


>  s  (x  ) 
iTi  3  i 


ll0,V'xilNlYxil! 


for  all  i=l,2....,M  and  j=l,2....,N. 


(14) 


The  code  book  resulting  from  such  a  cluster  design  may 
contain  a  high  degree  of  complexity  in  coder  implementation. 
This  arises  from  the  fact  that  an  exhaustive  search  of  the 
entire  code  book  must  be  made  for  the  nearest  output  vector 
to  the  input  vector  since  the  output  vectors  have  no  natural 
structure.  However,  at  the  price  of  suboptimality  and 
increased  code  book  storage  requirements,  one  can  induce 
certain  structures  upon  the  code  book  to  greatly  reduce  the 
search  time. 

Buzo,  Gray,  Gray,  Markel  [8]  use  a  binary  tree- 


structured  code  book  in  a  speech  coding  problem.  Juang, 
Wong,  Gray  [22]  show  a  more  general  M-ary  tree  structured 
code  book  design  offers  better  tradeoff  between  performance 


E.  LATTICE 


Another  structured  code  book  design  is  based  upon 
lattice  theory.  Gersho  [13]  thoroughly  describes  lattice 
quantizers.  The  quantizer  uses  a  set  of  output  points  which 
lie  in  a  bounded  region  of  a  lattice.  A  K-dimensional 
lattice  being  defined  by  any  nonsingular  KxK  matrix  (U)  so 
that  for  a  K-dimensional  column  vector  (m) ,  the  lattice  <A> 
is  given  by  the  set  of  all  vectors  of  the  form: 


=  U  m  (15) 

Notice  that  the  columns  of  U  are  lattice  points  and  form  a 
basis.  So  all  other  points  may  be  found  by  taking  linear 
combinations  of  the  basis  vectors  with  integer-valued 
coefficients.  An  example  of  a  two-dimensional  hexagonal 
lattice  quantizer  is  depicted  in  figure  11.  The  hexagon 
shapes  define  the  partition  regions  of  the  decision  space, 
and  the  points  correspond  to  the  representative  values  for 
each  partition  set.  Note  that  the  points  are  in  fact  the 
centroid  of  each  hexagon. 

Since  the  lattice  points  form  a  regularly  spaced  array 
of  points  in  K-dimensional  space,  a  lattice  quantizer  is  a 
uniform  quantizer.  Gersho  [19]  uses  this  fact  in  designing 
multidimensional  compandors,  much  as  those  of  the  one-dimen¬ 
sional  class  are  designed. 


intizer 


Intra-frame  DPCM  using  adaptive  quantization  has 
recently  been  studied  by  Zetterberg  [38]  for  two  types  of 
adaptive  quantizers  discussed  earlier  —  the  Jayant 
quantizer  and  a  forward  adaptive  quantizer  using  a  variance 
estimator.  Comparisons  were  made  and  tabulated  by  judging 
picture  quality  via  SNR  and  also  subjective  examination. 

They  observed  that  nonmoving  areas  should  use  a  fine 
quantizer  to  limit  granular  noise  to  low  levels  whereas 
moving  areas  need  a  large  dynamic  range  or  coarser  quantizer 
to  avoid  overload  by  large  prediction  values.  Entropy 
coding  was  then  used  in  order  to  meet  a  1  bit/pel 
requirement  of  European  standard  transmission  rate.  The 
basic  block  diagram  is  figure  12. 
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CHAPTER  THREE 


THE  ASPC  SYSTEM 

A  number  of  well  known  techniques  in  communcations 
theory  exist  for  compressing  data  from  a  source  signal.  Two 
methods  which  have  received  much  attention  are  the 
Dif feriential  Pulse  Coded  Modulation  (DPCM)  method  and  the 
transformation  method.  Both  exhibit  the  desired  data 
compression  performance,  yet  their  operation  is  quite 
different. 


I.  ffPCM  SYSTEMS 

DPCM  systems  operate  under  the  assumption  that  a 
degree  of  correlation  exists  among  a  stream  of  data  that  is 
to  be  transmitted.  This  assumption  is  the  basis  for  a 
scheme  whereby  a  prediction  of  each  datum  can  be  made  based 
on  the  previous  data.  In  other  words,  a  prediction  at  time 
k,  is  made  on  the  next  piece  of  data  to  be  transmitted,  at 
time  k+1,  based  upon  a  linear  combination  of  the  present  and 
previous  data.  So  that, 

rv 

S (k+1 I k, k-1 , . . . , k-n) =  \  a  S(k-i)  (16) 

i 

i=o 

where  S(k)  represents  the  data  stream  to  be  transmitted, 
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S(k+1  k, k-1 , . . . , k-n)  the  predicted  value  for  time  k+1  based 


on  the  previous  values  back  through  time  k-n,  and  (a  }  is  a 

i 

set  on  n+1  coefficients. 

This  prediction  is  then  subtracted  from  the  actual 
value  of  S  at  time  k+1,  S(k+1) ,  yielding  an  error  or 

"residual"  from  the  original  sequence,  called  E(k+1) : 


E(k+1)  *  S(k+1)  -  S(k+1  k, k-1, . . . ,k-n) 


This  residual  is  then  quantized  and  transmitted.  See  figure 
13-A.  If  the  predictor  and  quantizer  are  adequately 
designed  to  the  signals,  then  the  desired  data  compression 
is  achieved. 

The  receiver  stage  is  assumed  to  have  the  same 
predictor  as  the  transmitter,  but  since  it  will  only  receive 
a  quantized  estimate  of  the  true  residual,  the  transmitter 
must  be  restricted  to  operate  with  the  quantized  residuals 
in  generating  each  update  in  order  to  achive  parity  between 
the  two.  Therefore  the  quantized  residual  is  added  to  the 
prediction  to  obtain  the  source  estimate  to  be  used  in  later 
predictions.  So  in  fact, 


i(k+l  Ik,  k-1, . . .  ,k-n)  =  S  a  S(k-i) 


where 


S(k)  =  S(k |k-l, . . . ,k-n-l)  +E  (k)  (19) 

q 

and  the  quantized  residual  E  (k)  is 

q 

E  (k)  =  S(k)  -  S(k Ik-1, ... ,k-n-l)  +  Q  (k)  (20) 

q  n 

Thus  the  observed  or  estimated  sequence  is: 

S (k)  =  S(k |k-l, . . . ,k-n-l)  +  (21) 

S(k)  -  S(k Ik-1, . . . ,k-n-l)  +  Q  (k) 

n 

=  S(k)  +  Q  (k) 
n 

is  the  actual  sequence  S(k)  corrupted  by  quantization  noise, 

Q  (k)  . 
n 

At  the  receiver,  this  process  is  mirrored  by 
reconstructing  the  estimated  sequence  of  the  source  by 
adding  the  prediction  to  the  received  quantized  residual. 
See  figure  13-B.  Both  predictors  of  the  transmitter  and 
receiver  incorporate  an  n-cell  memory  register  for  the 
previous  n  values. 

So  if  one  neglects  the  effects  of  channel  errors  in 
the  transmission  and  any  external  noise  injected  into  the 
transmitter  or  receiver,  then  the  quantization  noise  Q 
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DPCM  a)  Transmitter  Loop 
b)  Receiver  Loop 


Figure  13 


proves  to  be  the  detrimonious  aspect  of  the  system  for  a 
particular  predictor  scheme.  Thus  for  a  fixed  prediction 
scheme,  the  quantizer  yields  the  limiting  influence  to  the 
system's  performance.  This  is  why  so  much  energy  has  been 
devoted  toward  the  study  of  quantizers  for  this 
communications  configuration. 

II.  TRANSFORMATION  SYSTEMS 

According  to  Pratt  [22] ,  the  notion  of  coding  and 
transmitting  the  Fourier  transform  of  a  monochromatic  image 
in  place  of  the  original  image,  was  introduced  by  Andrews 
and  Pratt  in  1968.  Images  of  this  type  tend  to  have  a  high 
degree  of  correlation  among  neighboring  pixels.  This 
correlation  causes  the  energy  distribution  in  a  two- 
dimensional  transform  of  a  monochromatic  image  to  be 
clustered  into  small  portions  of  samples.  It  was  discovered 
that  to  achieve  a  reduction  in  data  rate,  and  thus  bandwidth 
requirements,  only  the  high  order  transform  coefficients 
need  to  be  quantized  accurately  to  maintain  a  satisfactory 
reconstruction  of  the  image. 

The  normal  procedure  in  two-dimensional  image  coding 
is  to  apply  a  unitary  transformation  to  an  image  which  has 
been  divided  up  into  equal  size  square  blocks.  The  highest 
energy  transform  coefficients  in  each  block  are  finely 
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quantized  while  the  lower  energy  coefficients  are  grossly 
quantized.  These  quantized  values  are  then  encoded  and 
transmitted  block-by-block  to  the  receiver  which  performs 
the  appropriate  decoding  and  inverse  tranf ormation  to 
reconstruct  each  original  block,  and  hence  the  orginal 
picture.  See  figure  14. 

Since  this  concept  was  first  introduced,  many  have 
studied  the  effects  of  various  unitary  transformations  and 
block  sizes  on  the  image  reconstruction  quality  and 
bandwidth  reductions.  Some  of  the  transformations  examined 
include  the  Karhunen-Loeve,  Hadamard,  Haar,  Slant,  and 
Cosine  transforms.  The  Karhunen-Loeve  transformation, 
consisting  of  a  matrix  of  eigenvectors  of  the  data  to  be 
transformed,  provides  a  total  decorrelation  of  the  block  to 
give  the  lowest  distortion  of  any  of  the  transforms.  But 
the  problem  with  this  transformation  is  the  requirement  of 
an  exact  statistical  knowledge  of  the  image  block  to  be 
transformed,  which  is  normally  unavailable.  Even  if  it  is 
available,  the  computation  of  the  transform  is  slow  and 
complex.  Sub-optimal  approximations  of  the  Karhunen-Loeve 
transform  have  been  made  by  assuming  the  original  image  to 
be  of  a  known  specific  statistical  nature,  such  as  a  first- 
order  Markov  process.  If  this  type  of  assumption  is  made, 
then  the  computation  of  the  eigenvector  transformation 
matrix  becomes  fast  and  simple.  Unfortunately  this  is  much 
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inferior  since  the  original  images  would  not  usually  all  be 
of  this  nature,  nor  even  throughout  any  one  image. 

More  detailed  discussions  of  the  transformations  may 
be  obtained  in  Pratt  122] .  Of  the  remaining  transformations 
studied,  it  appears  that  the  Cosine  transform  gives  the  next 
best  performance,  irregardless  of  the  block  size,  followed 
by  the  Slant,  Hadamard,  and  Haar  transforms.  Algorithms  of 
these  transforms  exist  and  most  are  relatively  fast. 

III.  iiXBfilfl  SYSTEMS 

Both  the  DPCM  and  transformation  coding  systems  have 
their  advantages  and  disadvantages.  The  DPCM  system  is 
relatively  easy  to  implement,  but  produces  an  inferior 
reconstruction  as  compared  to  that  of  transformation  coding 
systems.  Unfortunately  transformation  coding  systems 
involve  a  more  complex  design.  By  combining  the  two  in  a 
particular  fashion,  it  has  been  discovered  that  the 
advantages  of  both  may  be  realized  while  minimizing  the 
disadvantages  of  each  system  alone. 

Habibi  [19]  proposed  that  a  system  be  comprised  of 
the  two  coding  schemes  in  such  a  manner  that  a  one¬ 
dimensional,  instead  of  a  two-dimensional,  transform  be 
applied  to  each  image  block  to  decorrelated  a  sequence  of 
the  picture  data  in  one  direction  throughout  the  block.  The 
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transform  coefficients  are  then  taken  across  the  block  and 
fed  into  the  input  of  a  DPCM  system  operating  in  the 
orthogonal  direction  to  that  of  the  transform.  One 
dimensional  transforms  are  fast  and  very  easy  to  implement 
and  combine  with  the  DPCM  system  to  give  quick,  reliable 
performance.  See  figure  15. 

As  briefly  mentioned  in  the  previous  chapter,  several 
variations  of  the  original  DPCM  concept  exist,  each 
providing  certain  advantages  in  particular  applications. 
One  of  the  most  useful,  especially  in  digital  image  coding, 
is  that  of  ADPCM  (Adaptive  Diffeiential  Pulse  Coded 
Modulation) .  With  only  a  small  increase  in  system 
complexity,  the  prediction  algorithm  of  the  usual  DPCM  loop 
can  be  made  to  adapt  to  the  statistical  variations  of  the 
source,  producing  more  accurate  predictions  for  new  data  and 
therefore  better  decorrelation. 

The  basic  idea  is  to  perform  an  analysis  of  the  data 
to  be  predicted  in  order  to  arrive  at  an  optimal  set  of 
coefficients  for  the  autoregressive  predictor.  This 
analysis  and  subsequent  coefficient  set  are  achieved  by 
introducing  a  learning  period  into  the  system.  The  data  is 
buffered  from  the  ADPCM  component  for  one  learning  period, 
causing  the  system  to  no  longer  be  real-time,  but  offset  by 
one  learning  period.  During  this  time  the  data  is  examined 
by  means  of  such  statistical  tools  as  a  Kalman  filter,  ARMA 


(AutoRegressive  /  Moving  Average)  model,  or  Stochastic 
Approximation  algorithm. 

For  a  detailed,  and  thorough,  analysis  and  comparison 
of  the  aforementioned  three  predictors  and  predictors  in 
general,  see  the  work  of  my  colleague  and  co- investigator  of 
Hybrid  /  DPCM  systems,  Stratigakis  [35] .  The  Kalman  filter 
is  well  known  from  Control  and  Estimation  Theory,  and  its 
derivation  is  prolific  in  the  engineering  literature  [6] , 
[26],  [27].  The  Kalman  filter  is  the  best  linear  minimum 
variance  estimator  and  for  the  case  of  gaussian  distributed 
signal,  noise,  and  initial  state,  it  is  the  best  of  all 
possible  linear  and  nonlinear  estimators  [27] . 

The  ARMA  model  is  a  method  of  determining  a  precise 
statistical  model  of  the  sequence  being  analyzed,  and  the 
estimates  of  the  autoregressive  coefficients  is  obtained  by 
solving  a  set  of  equations  involving  n  linear  combinations 
of  the  estimated  autocorrelations  of  the  sequence,  known  as 
the  Yule-Walker  equations  [51.  The  Stochastic  Approximation 
algorithm  is  similar  to  the  Kalman  filter  in  form  and 
operation.  It  is  extremely  fast  but  produces  an  inferior 
estimation. 

Unfortunately  these  algorithms  do  not  always  generate 
a  stable  set  of  coefficients  for  a  particular  line.  But  by 
making  use  of  the  partial  autocorrelations  of  the 
coefficient  set,  the  stability  of  the  shaping  filter  can  be 
determined  and  the  coefficient  set  adjusted  if  needed.  This 
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process  is  called  PARCOR  stabilization, 
and  adjustment  occur  as  follows: 


The  stabilization 


a)  The  partial  autocorrelations  of  the 

coefficients  (p  ,  i=l,n)  are  calculated, 
i 

b)  If  p  >  1  for  any  i,  then  those  coefficient 

i 

values  are  scaled  to  some  magnitude  less 
than  one  and  the  partial  autocorrelations 
are  recalculated. 

c)  The  process  iterates  until  all  partial 
autocorrelations  achieve  the  requirement  of 
their  magnitude  restriction. 


The  procedure  is  therefore  to  take  a  row  of 
transformed  picture  data,  apply  the  coefficient  set 
identification  algorithm  (Kalman  filter,  ARMA,  or  Stochastic 
Approximation) ,  and  stabilize  any  instabilities  in  the 
system  by  re-adjusting  the  coefficients  to  the  predictor  in 
the  ADPCM  loop.  Since  the  transmitter  will  make  its 
predictions  on  this  set  of  coefficients,  the  set  must  also 
be  made  available  to  the  receive’*  for  the  ADPCM  decoding. 
It  is  therefore  necessary  to  quantize  and  encode  this 
supplemental  information  into  the  data  stream.  Discussion 
of  this  quantizer  and  all  other  quantizers  for  the  hybrid 
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CHAPTER  FOUR 


Mx  -  Lloyd  Quantizer 


Optimization  of  a  quantized  representation  of  signals 
has  produced  considerable  interest  in  quantizer  research  and 
areas  where  rate-distortion  considerations  are  of 
importance.  Fundamental  to  most  of  the  research  to  date  on 
this  problem  are  the  works  of  Max  [25]  and  Lloyd  [24]  in' 
their  developments  of  digital  transmission  systems  which 
minimize  quantization  error. 

The  signal  to  be  quantized  is  assumed  to  be  an 
ergodic  process  of  known  probability  distribution.  By  con¬ 
sidering  the  input  signal,  S  ,  as  composed  of  an  indexed 

in 

class  of  sets  {R  ,  . ..,  R  },  where: 

1  N 


R  =  { x  x  <x<x  ,  x  =-oo,x  =  o°  }  (22) 

i  '  i  i+1  1  N+l 


to  which  a  corresponding  set  of  representative  values 

{y  .  ...,  y  }  is  assigned,  an  index  function  /(x),  -*«<x«», 
1  N 

on  the  partition  set  {R  .  ...,  R  }  is  defined  as: 

1  N 


yix)  =  [i  |  x<R  ,  l<i<N]  (23) 

i 
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and  a  sequence  o£  labels: 


,<t  )  =  J(S  (t  ) ) 
i  in  i 


representing  the  output  signal,  S  , 

out 

transmission. 


is  generated  for 


If  a  distortion  measure  on  the  quantization  mapping 
of  partition  sets  to  representative  sets  is  defined  as  the 
expected  value  of  some  differentiable  function  f(e),  where 
€  is  the  quantization  error,  and  the  known  probability 
density  of  the  input  signal,  S  ,  is  p(x) ,  then: 


D  »  E  [  f(S  -  S  )  ] 
in  out 

N  r 

=  V  /  f  (x  -  y, )  p ( x )  dx 


if 

1*1 4 


To  minimize  the  distortion  for  a  particular  fixed 

value  of  N,  necessary  conditions  may  be  obtained  by 

differentiating  the  distortion  (D)  with  respect  to  the 

partition  (R  )  and  representative  (y  )  sets  and  setting 
i  i 

<£■  equal  to  zero. 
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The  choice  of  f  is  usually  a  metric  function  so  that: 


f(0)  =  0 
f  ( x )  =  f(-x) 


If  it  is  required  that  Q(x)  be  monotonically  increasing, 
then  (28)  implies: 


x  -  y  =  x  -  y  i  =  2 ,  . . . ,  N 

i  i-1  i  i 


(yi  +  *l-!> 


i  =  2,  . , . ,N 


Usually  the  metric  function  is  chosen  to  be  the 


squared  metric: 


f  ( x)  =  x 


This  results  in: 


(y  +  y  ) 

i  i-1 


i=2 ,  . . . ,  N  ( 34a) 


y  =  2  x  -  y 
i  i  i-1 


i  =  2  ,  . .  .  ,  N 


(34b) 
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and 


i+1 

(x  -  y  )  p(x)  dx  = 
i 


(35) 


i  =  1 ,  . . . ,  N 


between  x  and  x 

i  i+1 

An  iterative  method  may  easily  be  implement 

solve  for  the  optimum*  in  the  mean— squared— error  sense* 

threshold  levels  {x  }  and  quantization  levels  {y  }  for  the 

i  i 

known  probability  density  function,  p(x) ,  via  equations  (34) 

and  (35)  .  Figures  of  the  optimum  thresholds  and 

quantization  levels  for  uniform  and  gaussian  probability 

distributions  appear  in  figure  16. 

Once  the  set  of  thresholds  and  quantization  levels 
have  been  determined,  any  sample  may  be  quantized  by 
determining  the  partition  set  to  which  it  belongs  and  the 
representative  value  for  that  set  by  equation  (23) . 

The  Max— Lloyd  quantizer  can  now  be  made  adaptive  in 
the  same  fashion  as  any  other  quantizer.  A  step-size 
adjustment  A (k)  is  defined  such  that  the  thresholds  and 
levels  of  the  quantizer  may  easily  be  scaled  to  some 
statistical  variation  estimate.  The  use  of  Max-Lloyd 


P(X) 

.V 

to 

* 

\\ 

\ 


.'*1 


*  \ 

“N 
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for  the  Gausslenly  Distributed  Cose 


quantizers  in  an  Adaptive  Stochastic  Picture  Coding  (ASPC) 
system  for  various  system  parameters  and  data  quantization 
has  been  investigated  and  a  detailed  discussion  of  this 
system  and  its  performance  will  appear  in  Chapters  Six  and 
Seven. 
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CHAPTER  FIVE 


As  previously  mentioned,  quantizing  blocks  or  vectors 

of  samples  in  lieu  of  scalar  quantization  approaches  a  more 

accurate  representation  of  the  source,  in  a  rate  distortion 

sense,  subject  only  to  vector  length.  By  defining  vector 
# 

quantization  as  a  surjective  mapping  of  K-dimensional  input 
vectors  to  corresponding  representative  vectors, 

K  K 

Q:R - ►  Y,  YCR  (36) 


one  can  get  an  intuitive  feel  for  the  concept  of  vector 

quantization.  A  typical  two-dimensional  vector  quantizer 

might  appear  as  in  figure  17.  Each  partition  represents  a 

decision  rule  in  the  vector  space.  These  partition  sets  and 

corresponding  representative  vectors,  X  ,  define  the  vector 

i  £ h 

code  book  where  each  X  is  the  centroid  of  the  i  partition 

i 

set. 

Since  this  K-dimensional  space  must  be  partitioned 
into  a  finite  number  of  regions,  it  is  appropriate  to 


examine  the  optimality  criteria  to  be  achieved  to  determine 
the  maximal  rate  distortion  performance.  Linde,  Buzo,  and 
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Gray  [23J  have  given  a  fundamental  paper  on  the  design  of 


vector  quantizer  algorithms.  Their  derivation  is  presented 

here  since  the  vector  quantizer  simulation  done  for  this 

report  was,  based  upon  one  such  algorithm.  If  an  N-level 

quantizer  is  defined  as  optimal  (or  globally  optimal)  if  it 

★ 

minimizes  the  expected  distortion,  then  x  is  optimal  if 

*  Q 

for  all  N-level  quantizers  x  : 

Q 


DU  )  «:  du  )  (37) 

Q  Q 


If  DU  )  is  only  a  local  minimum,  then  x  is  said  to  be 

Q  Q 

locally  optimum. 

Likening  the  distortion  measure  to  that  of  Max  and 
Lloyd,  the  measure  chosen  for  this  report  was  the  squared- 
error-  distortion: 


d  U,.x) 


(38) 


between  x,  the  sample  vector,  and  x ,  the  representation 

vector.  However,  as  with  Max  and  Lloyd,  the  derivation  is 

independent  of  the  particular  distortion  measure  chosen. 

Given  an  N-level  vector  quantizer  x  defined  by  the 

Q 

geometric  partitioning  set  P  =  {  P  ,  i=l,...,N}  and 

i 
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constrained  representation  alphabet  A  =  {  y  ;  i=l,...N}  then 

i 

the  expected  value  of  the  distortion  becomes: 


(39) 


where  E  [  d(i,y  )  Ixe  P  ]  is  the  conditional  expected 

i  '  i 
distortion  given  that  ieP  . 

i  /n 

If  a  particular  representation  alphabet  A  is  given, 

A 

but  a  partition  is  not,  then  a  partition  for  this  alphabet  A 

may  be  found  by  obtaining  the  mapping  of  each  a  into  the  ^ 

A  i 

A  minimizing  the  distortion  between  the  two  vectors, 

d(*,¥  ) •  So  the  desired  partition  is  found  by  choosing  each 
i 

vector  such  that  this  distortion  is  a  minimum  —  this  type 
of  clustering  is  well  known  to  those  familiar  with  pattern 
recognition  techniques  as  "nearest  neighbor"  clustering. 
Anytime  the  lowest  distortion  is  equal  among  a  number  of 
vectors,  any  convenient  tie-breaking  choice  may  be  made. 
The  partition  found  in  this  manner  yields: 


or 


* 

xeP  <====>  d(*,y  )  <  dU,y  j 

i  i  j 


D  (  {A  ,  P  (A) }  )  -  E  [  min  d  U,y)  ] 


(40) 

(41) 
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implying  that  when  compared  to  any  other  partition,  P: 

* 

D  (  {A  ,  P  (A)}  )  $  D  (  {£  ,  P(A>}  )  (42) 

* 

So  an  optimal  partition  set  P  may  be  found  for  any  fixed 

A 

representation  alphabet  A. 

Likewise  one  can  determine  a  best  possible 

* 

representation  alphabet  A  for  any  given  fixed  partition  set 
P.  If  the  distortion  measure  and  distribution  are  such  that 
each  set  P  with  nonzero  probability  in  K-dimensional 
Euclidean  space  contains  a  minimum  distortion  vector  £(P)  so 
that: 


E  [  dU,*(P))  |i<  P  ]  (43) 

=  min  E  [  d(*,y)  I  *  «  P  1 
II  1 


and  is  termed  the  centroid  of  the  set  P.  Therefore  no 

A 

representation  alphabet  A  =  {  y  ;  i=l,...,N  }  yields 

i 

smaller  distortion  than  the  alphabet, 


*(P)  =  (x(P  );  i=l , . . . ,N}  (44) 

i 


where  each  Jt(P.)  is  the  corresponding  centroid  of  partition 

P  since: 
i 
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N 

D  ({A  ,  P})  =  ]>  E  l  dU,y  )|i€p.]?r  <*tp  ,  (45) 

i=l 


N 


^  \  min  E  [  d(j£,U)  x^P  ]  Pr  (j£  e  P  ) 

^  i  i 


i=l  II 


=  D  ({£(P)  ,  P}) 


* 

Therefore  the  optimal  representation  alphabet  A  =  £(P)  for 
any  fixed  partition  set  P. 

By  using  a  combination  of  equations  (41)  and  (45) ,  a 

method  can  be  implemented  for  deriving  a  good  quantizer 

through  successive  iterations  of  these  constraints  given  any 

starting  quantizer  algorithms  using  this  idea  have  been 

proposed  for  both  known  and  unknown  distributions  of  vectors 

[23].  Given  an  initial  training  sequence  of  n  vectors,  an 

R 

M-level  vector  quantizer  with  M=2  ,  R=0,1,...  is  derived 

until  an  initial  guess  for  an  N-level  quantizer  is  obtained. 
The  optimal  N-level  quantizer  is  then  found  using  this 

guess.  The  algorithm  chosen  for  this  report  work  is  as 

\ 

fOllOWS: 
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(1)  Initialization:  Set  M=1  and  define 

A  (1)=5(A)  where  5(A)  is  the  centroid  of 
0 

the  training  set. 

(2)  Given  the  representation  alphabet  A  (M) 

0 

containing  M  vectors,  {y  ,i=l,...,M},  each 

i 

codevector  y  is  used  to  produce  two  code 
i 

vectors  y  +  /3and  y  -  /3  with  Q  a  fixed 
i  i 

perturbation  vector.  A  new  alphabet  A  of 

{y  +  /3  ,y  -/3,  i=l,...,M}  with  2M  vectors  is 

i  i 

now  obtained.  Let  M=2M. 

(3)  If  M=N,  then  the  operation  is  complete  and 

a  a/ 

the  alphabet  A  =A(N)  is  used  as  the  initial 
0 

reproduction  alphabet  for  the  N-level 
vector  quantizer  and  processing  continues 
with  (9).  If  M  is  not  yet  equal  to  N.  then 
proceed  to  (4)  . 

(4)  Set  ]»0  and  D  -  oo  and  eps>0  some 

~1  A  >*-> 

threshold.  Initialize  A  =A(M) . 

A  0 

(5)  Given  A  ={y  ,i=l,...,M},  find  the  minimum 

1  i 

distortion  partition  P(A  )={P  ;  i=l,...,N} 

1  i 

of  the  training  sequence: 


A  «  P  if  d^  ,y  )  <dU  ,y  ) 
j  i  j  i  j  k 

k 
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(6)  Find  the  average  distortion  for  step  1  by: 


D  =D  ( { A  ,  P  ( A  )  } ) 
1  1  1 


n-1 


2mln 

j=0 


j 


n 


rf  <D  “D  5/D  ^  eps  then  A  is  the 

1-111  1 
reproduction  alphabet  for  the  M-level 

quantizer .  Set  A  (M)*=A  and  continue  at 

0  1 

(2) .  if  the  (D  -D  )/d  >  eps  continue  to 

1-11  1 

(8)  . 

(8)  Find  the  optimal  representation  alphabet: 

£(P(A  ) )={£<P  ;  i=l,...,M} 

1  1 

for  P (A  ) 

1 

where  each  i(P  )  is  the  centroid  of 
partition  P  .  Set  A  *£(P(A  )).  Replace 

i  1+1  1 

1=1+1  and  continue  at  (5). 

(9)  Now  that  the  initial  N-level  vector 

quantizer  has  been  determined,  A  (N) ,  the 

0 

same  procedure  as  outlined  in  (5)  through 
(8)  is  performed  for  an  N-level  quantizer 
until  the  change  in  distortion  between 

code  book  calculations  drops  below  a 
specific  threshold  or  a  maximum  number  of 
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iterations  occurs.  The  resultant  code  book 
should  yield  a  performance  within  the 
designated  distortion  performance. 


The  scheme  in  which  this  algorithm  is  used  will  be  described 
in  the  following  chapter  on  system  configurations. 


CHAPTER  SIX 


SYSTEM  CONFIGURATIONS  AND  PROGRAM  DESCRIPTIONS 

I.  ASPC  /  Forward  Adaptive  Scalar  Max 

If  a  good  data  rate  versus  distortion  performance  is 
to  be  achieved  in  a  hybrid/DPCM  system  (also  known  as 
Adaptive  Hybrid  Picture  Coding  -  AHPC  -  or  Adaptive 
Stochastic  Picture  Coding  -  ASPC)  it  is  necessary  to  examine 
the  transmission  requirements  for  the  proper  quantizer 
selections.  A  block  diagram  of  the  proposed  transmitter  and 
receiver  is  shown  in  figures  18  and  19. 

This  is  the  basic  design  for  an  AHPC/ASPC  system  with 
an  autoregressive  predictor  (P)  and  a  forward  adaptive 
residual  Max  quantizer  (Q5) .  The  operation  is  as  follows: 

1)  The  image  desired  to  be  transmitted  is  broken 
into  strips,  each  of  size  BLK  X  NS  where  BLK  is 
the  desired  transform  block  size  and  NS  is  the 
number  of  samples  in  an  image  line. 

2)  A  series  of  NS  columns  are  passed  through  a 
forward  unitary  transformation  for  each  strip, 
producing  a  strip  of  transformed  picture  data 


Transmitter  with  Forward  Adaptive  Max  Quanti 


AS PC  Receiver  with  Forward  Adaptive  Max  Quantizer 


which  has  been  decor related  column-wise. 

Each  transformed  image  line  is  taken  (row-wise) 

singly  and  put  into  the  first  buffer,  BUFFERl . 

This  data  line  is  taken  in  parallel  from  BUFFERl 

and  the  optimal  predictor  coefficient  set  {a} 

for  that  line  is  obtained  through  a  Stochastic 

Adaptation  method  (Kalman  filter,  ARMA, 

Stochastic  Approximation)  which  also  produces 

the  mean  estimate  of  the  line,  »  .  Both  the 

^ S 

coefficient  set  {a  },  i=l,...,n,  and  mean 

i 

estimate  are  quantized  via  Ql  and  Q2  for 

transmission,  {a}  and  £  ,  and  made  available 

Q  SQ 

to  other  parts  of  the  system. 

Since  the  predictor  in  the  AD PCM  loop  bases  its 
prediction  on  n  past  values  contained  in  an 
internal  shift  register,  this  register  must  be 
initialized  at  the  beginning  of  each  line  to 
produce  a  sensible  prediction  of  the  first  few 
values.  The  quantized  signal  mean  estimate  is 
used  for  this  purpose  at  both  transmitter  and 
receiver  stages. 

The  values  are  fed  serially  form  BUFFERl  after 
the  coefficients  have  been  identified  for  that 
line.  While  the  system  transmits  this  data 
line,  the  next  row  of  transformed  image  data  is 
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passed  to  BUFFER1 . 

As  each  value  is  taken  from  BUFFER1  it  is  both 

passed  to  another  buffer,  BUFFER2,  and  an  ADPCM 

loop  with  perfect  quantization.  This  loop  is 

necessary  in  any  forward  adaptive  quantization 

scheme  since  a  variance  or  step-size  estimate 

must  be  obtained  prior  to  quantization.  The 

mean  of  the  residuals  is  also  obtained  so  that 

the  quantizer  may  be  adjusted  to  perform  about 

that  mean  estimate.  Since  the  residuals 

produced  by  the  Kalman  filter  and  ARMA  models 

should  be  gaussian,  an  N-level  scalar  Max 

quantizer  is  used  to  quantize  the  residuals  in 

Q5  with  an  expected  qaussian  normal  probability 

distribution  of  mean  zero  and  variance  one.  The 

quantizer  thresholds  and  representative  levels 

are  all  known  for  this  optimum  quantizer  a 

priori  to  both  transmitter  and  receiver.  These 

thresholds  and  levels  are  scaled  for  each  line 

of  residuals  by  a  step  size  determined  as  the 

standard  deviation  estimate  (obtained  from  the 

2 

quantized  variance  estimate  o  ) .  The 

eQ 

receiver  will  need  the  residual  mean  and 


estimate  for  each  line,  so  only  the  quantized 

2 

estimates,  n  and  a  ,  are  used  by  the 


transmission  stage  and  encoded  in  the  data 


stream  for  the  receiver.  A  single  variance 
estimate  for  each  line  is  assumed  sufficient  for 
the  statistical  variation  description  of  that 
entire  line. 

8)  The  values  are  taken  from  BUFFER2  and  the  final 

ADPCM  loop,  the  residuals  of  the  prediction,  e, 

are  acquired  and  quantized,  e  ,  and  passed  to 

Q 

the  coder.  The  quantized  residuals,  e  ,  are 

Q 

also  added  to  the  prediction,  r,  so  that  signal 

estimates  corrupted  by  quantization  noise,  s  , 

Q 

are  obtained  for  the.  predictor 1 s  n-stage  past 
value  shift  register. 


e  =  s  -  r  -  **  (46) 

SQ 

e  =  Q  (  e  )  =  e  +  q 

Q  n 

A 

s  =  e  +  r  +  M 
Q  Q  SQ 

=  s  +  q 

n 


9)  All  pertinent  data  required  by  the  receiver  for 
proper  reconstruction  of  the  image  is  then 
encoded  and  passed  through  a  digital  channel  to 
the  receiver. 
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10)  At  the  receiver  stage,  the  information  is 
decoded  into  its  various  components.  The 
quantized  signal  mean  estimate,  V-  SQ,  is  used  to 
initialize  an  n-stage  shift  register  in  a 
predictor  here,  and  the  quantized  coefficient 
set  for  this  data  line  is  provided  to  the 
predictor . 

11)  Each  prediction  value,  r,  is  added  to  the 

received  quantized  residual,  e  ,  and  adjusted  to 

Q 

the  appropriate  mean  to  obtain  the  signal 
estimate  of  each  transformed  pixel. 

s  =  e  +  q+  M  -  s  +  q  (47) 

Q  n  SQ  n 

12)  These  values  are  used  to  rebuild  the  original 
transform-domain  strips. 

13)  The  strips  are  finally  passed  through  an  inverse 
column  transformation  and  pieced  together  to 
reconstruct  the  original  image. 

Besides  the  quantizer  used  for  the  residual 
sequence,  Q5 .  four  other  quantizers  were  needed  for  the 
various  system  parameters,  none  of  which  were  to  be 
adaptive.  The  quantizers  were  uniquely  designed  for  each 
parameter.  The  quantizer  for  the  coefficient  set,  Ql,  was 


designed  as  an  8-bit  uniform  quantizer  of  the  limits  of  +/- 
1.75.  These  values  as  were  all  quantizer  parameters  for 
quantizers  Q1-Q4  were  determined  initially  by  estimation  and 
finalized  through  experimentation. 

The  quantizer  for  the  signal  mean,  Q2 .  was  originally 
designed  as  a  9-bit  fixed  Max-gaussian  quantizer  which  used 
two  sets  of  parameters.  One  set  was  used  with  those  image 
rows  occuring  at  the  beginning  of  each  block  of  transformed 
values  since  such  rows  contained  the  highest  energy  -  and 
thus  more  widely  variant  -  transform  coefficients.  The 
values  used  for  this  set  were  128  X  BLOCKSIZE  for  the  mean 
estimate  and  a  variance  estimate  of  about  130,000.  The 
estimation  of  128  X  BLOCKSIZE  for  the  signal  mean  came  from 
the  median  of  the  256  possible  gray  levels  in  a 
monochromatic  image  (128)  summed  over  the  block.  The  second 
set  was  used  for  all  remaining  data  lines  of  the  block  and 
consisted  of  a  zero  mean  estimate  and  a  variance  estimate  of 

i 

350.  It  was  later  determined  that  for  a  particular 
transformation,  these  signal  mean  values  could  be 
approximated  to  eliminate  the  necessity  of  quantizing  and 
transmitting  this  parameter  (the  approximation  being  known  a 
priori  to  both  transmitter  and  receiver)  to  actually  improve 
the  system's  performance.  In  the  case  of  the  Cosine 
transformation,  signal  mean  estimates  of  2750  for  lines  at 
the  beginning  of  a  block  and  0  for  all  other  lines  were 
used.  This  elimination  of  the  need  to  transmit  the  signal 
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mean  gave  a  data  rate  savings  of  9  bits  per  line.  These 
values  would,  of  course  be  different  for  other 
transformations.  Exact  system  data  rates  will  be  calculated 
in  the  next  chapter. 

The  residual  mean  estimate  was  assumed  gaussian  over 
the  entire  image  and  thus  a  time  invariant  Max  quantizer, 
Q3,  was  used  here.  It  was  of  6-bits  in  resolution  with  mean 
of  zero  and  variance  of  150.  The  residual  variance  estimate 
quantizer,  Q4,  was  uniform  from  0  to  40,000  for  the  ARMA  and 
Kalman  filter  estimators,  and  0  to  100,000  for  the 
Stochastic  Approximation  estimator.  Since  this  parameter 
yielded  the  step-size  estimate  for  the  adaptive  residual 
quantizer,  a  resolution  of  9  bits  was  used.  The  reason  for 
the  wider  distribution  of  the  variance  for  the  Stochastic 
Approximation  routine  is  that  it  performs  only  a  quick 
estimation  of  the  predictor  coefficients,  thereby  causing 
the  generated  residual  sequence  to  be  of  a  more  variant 
nature. 

II.  ASPC  /  Backward  Adaptive  Scalar  Max 

As  described  before,  a  system  can  be  used  such  that 
no  extra  quantizer  information,  such  as  a  step-size 
adjustment,  is  needed.  This  system  would  use  the  backward 
adaptive  quantizer  which  acquires  its  step-size  adjustment 
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from  previously  quantized  values  readily  available  at  both 
transmitter  and  receiver  stages. 

It  was  seen  in  Chapter  Two  that  there  are  as  many 
forms  of  step-size  algorithms  for  the  backward  adaptive 
quantizer  as  for  the  forward  adaptive  quantizer.  The  one 
selected  for  this  study  was  that  of  equation  (4),  reprinted 
here  for  convenience: 

i-  “  - 

1  V7  i-1  . 

A  (k)  =  2Lal  eQ(k-i)  |  (48) 

a 

2 

This  particular  formula  for  the  step-size  was  chosen  for  its 

similarity  to  the  autoregressive  prediction  being  performed 

in  the  ADPCM  loop.  Choices  of  scalars  a  and  a  are 

1  2 

determined  experimentally  for  particular  system.  The 

simulation  performed  for  this  report  used  values  of  a  =0.7 

1 

and  a  =1.  Future  inquires  into  this  algorithm  could 

2 

possibly  discover  more  suitable  parameters,  but  these  were 
judged  significant  and  typical  for  this  comparison  of 
quantization  schemes. 

With  a  backward  adaptive  quantizer,  the  transmitter 
of  figure  18  is  modified  somewhat  since  a  prior  residual 
variance  is  no  longer  needed.  However  the  residual  mean  is 

i 

still  required  for  proper  orientation  of  the  residual 

i 

quantizer.  The  system  does  require  additional  circuitry 
since  an  estimate  must  be  made  from  previous  quantized 
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Transmitter  with  Backward  fidaptive  Max  Quanti 


values.  The  resulting  transmitter  block  diagram  is  depicted 
in  figure  20.  The  receiver  block  diagram  remains  as 
shown  in  figure  21. 

III.  ASPC  /  Vector  Quantization 

An  extension  of  the  system  using  scalar  quantization 
techniques  of  figures  18  and  19  results  in  the  vector  AHPC 
system,  called  VASPC,  as  depicted  in  figures  22  and  23. 
The  operation  of  the  VASPC  system  is  nearly  identical  with 
that  of  the  ASPC  system  —  with  the  obvious  exception  of  the 
vector  quantization. 

The  original  image  is  still  split  into  strips  of  the 
transformation  blocksize  by  the  image  width  and  forward 
transformed  columnwise.  But  now  instead  of  applying  ADPCM 
on  a  single  line  at  a  time*  the  system  has  the  ability  to 
handle  any  number  of  lines  at  once.  Since  the  usual  image 
column  length  is  a  multiple  of  2,  and  the  transformation 
blocklength  also  a  multiple  of  2,  say  16,  the  vector 
dimensions  examined  are  K  =  2,  4,  8,  16.  It  must  be  kept  in 
mind  that  any  practical  size  can  be  specified,  but  any 
additions  in  dimension  increases  the  corresponding  vector 
code  book  /  alphabet  as  well  as  system  complexity. 

The  coefficient  set  for  each  line  is  determined  in 
precisely  the  same  fashion  as  the  scalar  AHPC  case.  The 
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first  vector  ADPCM  loop  figured  in  the  diagram  serves 

another  purpose  than  that  of  the  scalar  case.  Recall  that 

previously  this  particular  loop  was  for  the  determination  of 

an  estimate  of  the  residual  variance  of  each  transformed 

image  line.  In  the  VAHPC  system  this  loop  must  aid  in 

determining  the  whole  vector  quantizer  code  book.  K  lines 

are  operated  on  at  once  in  parallel  and  are  considered  to  be 

th 

NS  (number  of  samples  per  line)  vectors  in  the  K 
dimension.  A  corresponding  predicted  K-dimensional  vector 
is  produced  for  each,  and  the  resultant  residual  vectors  are 
collected  to  form  a  training  set  for  the  design  of  the 
vector  quantizer. 

£  =  £  -  B  i=l,...,NS  (49) 

i  i  i 

£  =  t  £  }  (50) 

i  i=l , . . . , NS 

This  set  of  training  vectors  is  then  fed  into  the 
design  algorithm  to  produce  a  quantizer  for  this  strip  of  K 
X  NS  residuals.  Since  this  code  book  must  be  available  to 
the  receiver,  a  mean  estimate  vector  and  a  variance  estimate 
vector  for  this  code  book  are  generated,  and  a  K-dimensional 
Max  gaussian  quantizer  scaled  to  these  statistical  vectors 
is  employed  to  quantize  and  transmit  the  code  book  to  the 


receiver.  For  small  vector  dimensions  the  scheme  of  design 
and  transmission  could  yield  an  enormous  overhead  if  it  is 
done  for  every  NSF  (number  of  image  lines)  /  K  size  strip. 
Instead,  after  each  code  book  is  designed,  each  codevector 
is  compared  with  each  codevector  of  the  previous  K  X  NS 
strip  code  book.  A  distortion  measure  is  obtained  from  this 
comparison,  and  the  present  code  book  is  shuffled  such  that 
those  vectors  achieve  a  minimum  distortion.  Those  falling 
below  a  specific  threshold  are  not  transmitted.  No  vector 
code  book  is  necessarily  similar  to  the  previous  one,  and 
therefore  a  vector  which  need  not  be  transmitted  does  not 
necessarily  have  the  same  index  as  the  vector  it  resembles 
of  the  previous  code  book.  Because  of  this  fact,  it  is 
necessary  to  give  an  index  to  the  old  code  book  for  which  a 
vector  is  being  replaced.  This  indexing  itself  constitutes 
an  overhead,  so  a  decision  must  be  made  as  to  whether  it  is 
better  to  transmit  a  whole  new  code  book,  a  partial  code 
book,  or  retain  the  previous  code  book.  This  decision  also 
causes  a  small  2  bit  /  code  book  overhead,  but  is  unchanging 
throughout  the  picture.  This  decision  takes  the  form  of: 

B 

K  *  NBE  *  2  =  B  *  NUMVTR  +  K  *  NBE  *  NUMVTR  (51) 

=  (B  +  K*NBE)  NUMVTR 
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where : 


K  =  Vector  dimension 

NBE  =  Number  of  bits  used  to  quantize  each 
element  in  each  vector. 

NUMVTR  =  Number  of  vectors  needing  to  be 
transmitted  . 

So  for  a  4  bit  vector  quantizer,  B=4  (16  codevectors)  and 
the  decision  becomes: 


MU 

2 


Transmit  wkQ.l.e  book  jJLl 
NUMVTR  >12 


4  NUMVTR  >13 

8  NUMVTR  >14 

16  NUMVTR  >  15 


and  for  a  16-dimensional  quantizer,  a  partial  code  book  is 
transmitted  even  if  only  one  vector  is  similiar  to  a 
previous  vector.  If  the  whole  book  is  to  be  transmitted, 
the  indexing  is  dropped.  For  a  comprehensive  data  rate 
calculation,  refer  to  Chapter  Seven. 

IV.  ASPC  /  Forward  Adaptive  Scalar  Max  with  "Zeroing" 


The  last  configuration  considered  is  actually  a 
degenerative  form  of  the  first  ASPC  configuration  with  a 


forward  adaptive  scalar  Max  gaussian  residual  quantizer. 
This  system  was  simulated  only  to  give  an  approximation  of 
AHPC  performance  at  non-integer  bit  rates  (neglecting 
overhead) .  Since  there  is  no  way  to  design  a  Max  quantizer 
with  0.5.  1.5.  2.5  or  any  non-integer  bit  resolution,  a 

method  called  conditional-replenishment  or  "zeroing"  is  used 
to  get  a  more  generalized  picture  of  the  scalar  AHPC  rate 
distortion  performance.  Even  if  one  could  make  such  a  1.5 
bit  quantizer  how  would  one  transmit  1.5  bits? 

Since  only  a  performance  estimate  is  desired,  it 
is  suitable  to  design  for  an  effective  performance.  The 
procedure  is  simply  to  not  transmit  portions  of  the  residual 
sequence,  thus  reducing  the  data  rate.  So  a  residual  of 
zero,  or  a  perfect  prediction  is  assumed,  and  the  subsequent 
error  in  the  next  prediction  is  expected  to  be  corrected  in 
the  residual,  which  is  transmitted.  Residual  sequences 
might  then  appear: 

e  0  e  0  e  0  . . . 

13  5 

e  0  0  e  0  0  e  0... 

14  7 

e  000e  0  0  0  e  0... 

15  9 
etc. 
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The  effective  bit  rate  is  then: 


B  Log  N  (#  bits  in  res.  qntzer) 

2 


Number  of  Number  of  sharing  residuals 

sharing  res. 

Of  course  the  performance  is  quite  inferior  at  all 
data  rates  to  the  normal  AHPC  system  since  the  coefficient 
identification  routines  are  designed  for  a  different 
sequence.  Instead  of  the  n-stage  past  value  shift  register 
containing  samples  perturbed  only  by  quantization  noise: 

s  =  e  +  r  =  e+r  +  q  (52) 

Q  Q  n 

there  is  the  additional  burden  of  a  previous  prediction 
error : 

A/  ,  , 
s  =  e  +  r  +  e  =  e  +  r+  q  +  e  (53) 

Q  Q  n 

and  the  prediction  scheme  uses  coefficients  which  have  not 
been  optimized  for  this  sequence.  But  the  routine  still 
yields  a  general  "feel"  of  the  scalar  AHPC  rate  distortion 
performance  at  odd  data  rates. 


v.  Program  Descriptions 

The  computer  programs  used  for  the  simulations 
presented  here  are  slightly  lengthy  but  have  been  designed 
in  as  modular  a  fashion  as  possible  to  achieve  a  flexibility 
in  design  changes  and  alternatives. 

The  transformation  program  consists  of  a  driving 
program  and  a  set  of  transform  subroutines.  The  image  to  be 
used  is  selected  from  the  disk  on  which  it  has  been  stored, 
and  provided  with  the  desired  transformation,  blocksize,  and 
transformation  direction(s)  (row  or  column),  the  image  is 
transformed  and  the  resultant  transformed  image  is  placed  on 
disk.  All  disk  operations  are  performed  with  a  subroutine 
called  DISKIO  which  handles  the  image  files  a  line  at  a  time 
in  any  desired  manner. 

The  main  transmitter  /  receiver  program  consists  of  a 
small  main  driver  segment  to  read  the  original  transformed 
image  and  write  the  received  and  reconstituted  transformed 
image  to  disk.  The  main  subroutine  is  called  COMPMC  which 
performs  the  simulation  of  ADPCM  transmission,  reception, 
and  reconstruction.  Minor  changes  are  introduced  into  each 
run  of  COMPMC  which  performs  the  simulation  of  ADPCM 
transmission,  reception,  and  reconstruction  of  the 
transformed  image  data.  Minor  changes  are  introduced  into 
each  run  of  COMPMC  for  the  specific  coefficient  prediction 
schemes  and  all  scalar  quantizers. 


Subroutines  which  are  also  required  include: 

MAXQTZ  -  subroutine  for  designing  any  N-level 
scalar  Max  /  Lloyd  quantizer  using  the 
techniques  derived  by  Max. 

PARSBL  -  subroutine  for  the  partial  auto¬ 
correlation  stabilization  of  the  predictor 
coefficients. 

CENTRD  -  subroutine  for  finding  the  centroid  level 
in  a  probability  distribution  between  two 
thresholds  for  MAXQTZ. 

INTGRL  -  a  16-point  gaussian  quadrature  integration 
subroutine  required  by  both  CENTRD  and 
MAXQTZ . 

Functions  used  are: 

UNIFRM  -  function  to  provide  an  N-level  unifromly 
distributed  quantization  between  two 
specified  intervals. 

FUNC  -  function  to  generate  a  desired  probability 
denstity  function. 

SQERR  -  function  to  generate  the  squared  error 
distortion  of  a  value  in  a  particular 
probability  distribution. 

QNTZ  -  function  to  quantize  a  sample  via  Max's 

definition  given  a  variance  estimate  and 
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the  quantization  thresholds  and  levels 
produced  by  MAXQTZ .  The  step-size  scale 
factor  is  determined  as  the  square-root  of 
the  variance  estimate  (std  deviation) . 

The  VASPC  vector  program  uses  a  subroutine  called 
COMPV  quite  similar  to  COMPMC.  Special  provision  had  to  be 
made  for  the  vector  quantizer  but  the  general  flow  of  the 
routine  was  the  same.  Besides  the  aforementioned 

subroutines,  COMPV  requires: 

VQDSN  -  subroutine  for  random  vector  quantizer 
design  based  on  a  given  training  set. 

VECQNT  -  subroutine  to  determine  the  representative 
or  quantized  value  of  a  vector  via  the 
"nearest  neighbor"  approach. 

COMPMC  and  COMPV  both  require  subroutines  for  the 
desired  coefficient  identification  routine,  with  the 
exception  of  the  Kalman  filter  routine  which  had  been 
directly  incorporated  into  one  version  each  of  COMPMC  and 
COMPV.  Other  identification  routines  require: 

FT AUTO  -  IMSL  subroutine  for  finding  the  auto¬ 
correlations  of  a  given  sequence. 

FT ARPS  -  IMSL  subroutine  for  generating  a  set  of 
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coefficients  for  a  designated  Auto- 
Regressive  Moving  Average  (ARMA)  model 
given  the  autocorrelation  produced  by 
FT AUTO. 

STAPRX  -  subroutine  for  determining  the  prediction 
coefficients  via  the  Stochastic 
Approximation  algorithm. 

Finally,  the  received  and  reconstructed  transformed 
images  are  inverse  transformed  with  the  same  original 
transform  program  package  and  the  final  image  is  placed  on 
disk.  Programs  were  also  used  to  calculate  the  signal-to- 
noise  ratio  between  the  original  and  final  images  and 
generate  the  corresponding  error  histograms.  These  programs 
are  SNR  and  HISTl,  respectively.  General  purpose  routines 
were  used  for  the  transfer  of  images  form  disk  to  tape  and 
vice  versa  for  subjective  viewing  on  the  COMTAL  Vision  One 
Image  Processing  system  of  the  Visual  Interactive  Processing 
laboratory,  a  part  of  the  University  of  Arkansas' 
Engineering  Experiment  Station.  Flowcharts  and  program 
listings  are  included  in  Appendix  B,  and  the  reader  is 
referred  to  this  appendix  for  further  information  on  the 
details  of  the  simulations. 


CHAPTER  SEVEN 


RESULTS  AND 


A  number  of  orthogonal  transformations  exist  which 
could  have  been  used  in  the  simulations,  but  the  aim  of  this 
study  was  not  the  comparison  of  transforms  for  AHPC.  The 
simulations  mentioned  in  Chapter  Six  were  performed  for  two 
types  of  transforms:  the  Hadamard  and  the  Cosine.  These 
two  were  chosen  for  their  frequent  use  in  other  work  in  this 
area  of  research  and  widespread  use  in  current  scientific 
publications.  Further,  the  Hadamard  represents  a  transform 
of  average  complexity  and  performance,  while  the  Cosine's 
performance  appears  to  excel  in  transformation  coding 
systems  [22]  .  Finally  these  were  the  two  transforms  used  by 
Habibi  [19]  if  his  original  work  on  Hybrid  /  DPCM  systems. 
Only  the  Cosine  was  used  after  it  was  determined 
experimentally  to  perform  better  for  AHPC  systems  than  the 
Hadamard  for  any  coefficient  identification  scheme. 

Two  images  were  used  as  originals  in  the  initial 
stages  of  the  system  simulations.  These  are  :ie  128  X  128 
monochrome  image  of  the  well-known  MIT  girl  used  in  much  of 
the  literature,  and  the  256  X  256  monochrome  image  of  an 
aerial  photograph  of  a  city.  See  pictures  <1>  and  <2>  in 
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Appendix  A.  ASPC  simulations  were  performed  on  each  image 
with  both  transformations  and  with  and  without  PARCOR 
stabilization  to  determine  the  best  overall  coefficient 
identification  /  transform  scheme  to  use  in  the  quantizer 
comparison.  It  was  determined  that  the  best  overall 
performance  was  given  for  the  cosine  transformed  city  image 
with  a  set  of  3  predictor  coefficients,  controlled  w-ith 
PARCOR  stabilization.  See  figure  .24.  All  ASPC  simulations 
on  this  table  were  performed  with  a  transform  block  size  of 
16,  a  3-bit  forward  adaptive  residual  Max  gaussian 
quantizer,  an  8-bit  uniform  coefficient  quantizer,  a  6-bit 
Max  gaussian  residual  mean  estimate  quantizer,  a  9-bit 
uniform  residual  variance  estimate  quantizer,  and  a  9-bit 
Max  gaussian  signal  mean  quantizer.  These  yield  fairly  high 
data  rates,  but  a  good  objective  viewpoint  of  system 
performance. 

The  cosine- transformed  city  image  was  then  used  in 
determining  the  best  overall  predictor  performance.  As  seen 
from  figure  25,  the  Kalman  filter  and  ARMA  model  resulted 
in  almost  identical  rate  distortion  performance  and  were 
significantly  superior  to  the  Stochaotic  Approximation 
algorithm.  The  Kalman  filter  was  then  chosen  to  perform  the 
coefficient  identification  for  the  quantizer  comparisons. 
The  Kalman  filter  adaptive  identification  process  is 


summarized  [21]  : 


Message  Model  A(k+l)=A(k)  (53) 

T 

Observation  Model  s(k+l)=R  (k) A ( k ) +e(k+l) 

A  aP 

Identification  eqn  A(k+1) =A(k) +K (k+1) e (k+1) 


V„(k)R  (k) 

A  P 

Gain  equation  K (k+1)  =  - - - 

V  +R~  <k)V„(k)R  (k) 
Z  p  A  p 

Variance  eqn  V^(k+1) =[I-K(k+1)R  T(k)  ] 

A  p 

*  V„(k) 

A 

where  at  time  k: 

A ( k )  -  Coefficient  vector 

s(k)  -  sample  value 

e(k)  -  residual  value 

R  (k)  -  past  value  vector 
PT 

R  A(k)  -  predicted  value 
P 

K(k)  -  Kalman  gain 

V  -  error  variance  =  100 
Z 

V^Ck)  -  a  posteriori  error  variance 
A 


A  graph  depicting  the  performance  of  each 
quantizer  routine  is  shown  in  figure  26.  The  corresponding 
images  appear  in  pictures  <  3  -  27>.  As  predicted  by  the 
theory,  the  vector  quantizer  produces  the  best  overall 
subjective  and  quantitative  performance,  followed  closely  by 
the  forward  adaptive  scalar  version  and  a  quite  inferior 
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backward  adaptive  scalar  version. 

Although  there  does  not  appear  to  be  a  consensus  on 
the  calculation  of  a  signal-to-noise  ratio  performance  for 
images  which  would  greatly  aid  the  comparisons  and 
evaluations  of  new  and  different  techniques  —  a  ratio  of 
signal  variance  to  noise  variance  was  arrived  at  as  a  good 
critical  evaluation.  The  data  rate  calculations  are  as 
follow : 


Forward  Adaptive  Scalar  Max  Residual: 

Coeff.  reqs.  256  lines  *  3  coef/line  *  8  b/coef 

256  lines  *  9  b/line 
256  lines  *  6  b/line 
R  bit 8/ pel 


Res.  Var. 
Res.  Mean 
Res.  Qntzr 


==>  R  b/pel  +  0.15  b/pel  overhead 


Backward  Adaptive  Scalar  Max  Residual: 

Coeff.  reqs,  256  lines  *  3  coef/line  *  8  b/coeff 

256  lines  *  6  b/line 
R  bits/pel 


Res.  Mean 
Res.  Qntzr 


=->  R  b/pel  +  0.12  b/pel  overhead 
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Vector  Quantizer: 


No  Code  book  transmission 

2  *  Number  of  no  transmissions  (decision) 

Whole  Code  book  transmission 

2  *  Number  of  full  transmissions  (decision) 

Log  N  *  K  *  7b/elem  *  Number  full  trans. 

2 

K  *  9bits  *  Number  full  trans.  (variance) 

K  *  6bits  *  Number  full  trans.  (mean) 

Partial  Code  book  transmission 

2  *  Number  of  partial  transmissions  (decision) 

Log  N  *  K  *  7b/elem  *  Number  part,  trans. 

2 

K  *  9bits  *  Number  part,  trans.  (variance) 

K  *  6bits  *  Number  part,  trans.  (mean) 

==>  Actual  data  rate  must  be  calculated 
interactively  during  system  simulation. 

These  systems  show  the  feasibility  of  various  digital 
image  coding  and  transmission  schemes.  Vector  quantization, 
although  a  relatively  new  technique,  appears  to  have  great 
promise  in  any  such  data  reduction  system.  As  seen  here  the 
AHPC  system  can  yield  good  results  at  low  data  rates. 
Pictures  and  error  histograms  supporting  this  conclusion 
appear  in  Appendix  A.  All  signal-to-noise  ratio 


calculations  were  made  as  a  ratio  of  the  variances: 


Future  Research 

The  AH PC  system  has  been  shown  to  give  good  rate 
distortion  performance.  The  limiting  factor  has  also  been 
shown  to  be  the  quantizer  performance.  Optimization  of  any 
initial  parameters  such  as  those  used  in  the  identification 
routines  may  improve  performance,  as  well  as  determining  any 
better  step-size  estimators.  Any  future  research  endeavor 
into  AHPC-type  systems  should  therefore  include  any  new 
thoughts  as  to  the  possible  evolution  of  quantization 
theory.  Vector  quantization  research  appears  the  most 
lucrative  for  low  data  rate  performances.  Although  other 
forward  adaptive  and  backward  adaptive  quantizer  algorithms 
may  exist  which  yield  better  results  than  those  presented 
here,  these  are  figured  typical  and  any  such  scheme  appears 
limited  by  the  one-dimensional  scalar  quantization  in 
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adaptive  systems. 

One  possibility  of  future  research  efforts  is  the 
incorporation  of  the  AHPC  system  into  a  real-time 
transmission  of  sequential  images.  Special  quantization 
techniques  may  be  employed  in  such  systems  which  take 
advantage  of  previous  image  frame  information. 
Considerations  of  any  other  interframe  knowledge  could 
improve  an  intraframe  performance  of  AHPC. 

A  problem  which  might  also  be  scrutinized  is  that  of 
an  introduction  of  noise  into  the  system.  What  are  the 
effects  of  channel  errors  on  the  receiver's  reconstructed 
image?  How  does  initial  noise  in  the  original  image  and 
noise  in  the  image  in  the  transform  domain  affect  system 
reliability?  Perhaps  various  filtering  techniques  could  be 
applied  at  points  in  the  system  to  limit  image 
deterioration.  Only  more  research  can  provide  the  answers 


to  these  questions. 
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RETURN 
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***  Main  Driver  Program  for  Transforms  *** 


EXTERNAL  HADGEN,  COSINE ,  HAARGN,  KLGEN,  FOURGN,  DUMMY,  SLAGEN 
INTEGER  HAD/ 4HHAr  / , HAAR/4HHAAR/ , COS/4HCOS  / , OTYPE, OBSZ , 

1  SLAN/ 4HSLAN/ , KL/ 4HKL  / , FOUR/ 4HFOUR/ , ID/4HID  / 

COMMON/ CO DPRM/  IR, IC, IBSZ , NTYPE , OBSZ , OTYPE , IBLKSZ 
COMMON  /CORE/  TEST (1024) 

COMMON/FT/  KFFT 
REAL  XIN ( 256 , 256 ) 

COMMON/ B/  IDA, N 
NN=1 6 
ISIZE=16 
KFFT=0 
IBSZ=16 
OBSZ =16 
WRITE(6,200) 

200  FORMAT ( ' 0 ' ,//, ' A  AAAAAAAAAAAAAAAAAAAA') 
WRITE (6,201) 

201  FORMAT ( ' 0 ' ,//) 

DO  3000  KKK=1 , 1 

READ (5,10) ( ITYPE , M, N, IBLKSZ, IR, IC, NTYPE, LIPT, OTYPE, LOPT) 
10  FORMAT ( IX, A4, IX, 2 15, IX, 13 , IX, 2 12 , IX , A4 , 13 , IX, A4 , 13 , IX , II) 
IF (ITYPE. EQ. HAD)  GO  TO  100 
IF (ITYPE. EQ.HAAR)  GO  TO  110 
IF (ITYPE. EQ. COS)  GO  TO  120 
IF ( ITYPE. EQ. SLAN)  GO  TO  130 
IF ( ITYPE. EQ.KL)  GO  TO  140 
IF ( ITYPE. EQ. FOUR)  GO  TO  150 
IF (ITYPE. EQ. ID)  GO  TO  160 
20  WRITE (6, 25) ITYPE 

25  F0RMAT(1H1, //, 133 (1H?) ,//,6X,A4, IX, 'TRANSFORM  NOT  ALLOWED 
1  // ,133 (1H?) ) 

100  CALL  XFORMB (HADGEN, M,N, LIPT, LOPT, XIN) 

GO  TO  170 

110  CALL  XFORMB (HAARGN, M,N, LIPT, LOPT, XIN) 

GO  TO  170 

120  CALL  XFORMB (COSINE, M,N, LIPT, LOPT, XIN) 

GO  TO  170 

130  CALL  XFORMB ( SLAGEN, M,N, LIPT, LOPT, XIN) 

GO  TO  170 

140  GO  TO  20 

150  KFFT=1 

160  CALL  XFORMB (DUMMY, M,N, LIPT, LOPT, XIN) 

170  WRITE (6, 180) ITYPE, IR, IC 

180  F0RMAT(1H1,//,133(1H$) ,//,45X,A4, IX, 'TRANSFORM  COMPLETED' 

1  / , 45X, ' IR= ' ,12,/ ,45X, 'IC=' ,I2,//,133(1H$) ) 

888  REWIND  3 

IF(KKK.NE.l)  GO  TO  3000 
DO  72  1=1,256 

72  READ (3) (XIN(I,J) ,J=1,256) 

REWIND  3 
3000  CONTINUE 


***  Subroutine  for  Transformations  *** 


SUBROUTINE  XFORMB (RTN, M, N, LIPT, LOPT, DATA) 

Q  it  it  it  it  ir  it  it  it  it  it  it  *  it  it  it  *  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it 

C  PURPOSE:  MAIN  PROGRAM  /  TRANSFORM  SUBROUTINES 
C  INTERFACE. 

C  RTN  =  TRANSFORMATION  TO  BE  USED. 

C  MXN  =  IMAGE  DATA  SIZE  (M  BY  N) 

C  LIPT  =  DATA  LINE  INPUT  DISK  UNIT  NUMBER. 

C  LOPT  =  DATA  LINE  OUTPUT  DISK  UNIT  NUMBER. 

C  DATA  =  IMAGE  DATA 

Q  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it 

C  USER  NOTE:  BE  CERTAIN  TO  INCLUDE  ALL  NECESSARY  COMMON 
C  STATEMENTS  IN  MAIN  PROGRAM. 

£******************************************************** 
EXTERNAL  RTN 
REAL  A(1024) ,DATA(M,N) 

INTEGER  OBSZ , OTYPE 
COMMON/FT/  KFFT 

COMMON/ CODPRM/  IR, IC, IBSZ , NTYPE, OBSZ , OTYPE, IBLKSZ , I FLAG 
MAXRC  =  (1+KFFT) *IBLKSZ*N 
IF (MAXRC.GT. 32800)  GO  TO  500 
IF (KFFT)  10,10,20 

10  CALL  BLKFRM (RTN, M,N, LIPT, LOPT, DATA) 

RETURN 

20  CALL  CBKFRM (RTN, M,N, LIPT, LOPT, IBLKSZ, DATA) 

RETURN 

500  MXCORE  - (MAXRC  *  4)/1000 
WRITE (6. 501)  MXCORE 

501  F0RMAT(1H1,132(1H*) ,//,16X, 115, 'K  REQUIRED,  LINE  BY  LINE 
* ' NOT  YET  IMPLEMENTED' ,//, 133 (1H*) ) 

RETURN 


***  Subroutine  to  Break  Image  into  Blocks 


★  *  ★ 


SUBROUTINE  BLKFRM (RTN,  M, N, LIPT, LOPT, DATA) 

C**************************************** **************** 

C  PURPOSE:  IMAGE  BLOCK  FORMATION  AND  MANIPULATION. 

C  RTN  =  TRANSFORMATION  TO  BE  USED. 

C  MXN  =  IMAGE  SIZE  (M  BY  N) 

C  LIPT  =  LINE  INPUT  DISK  UNIT  NUMBER. 

C  LOPT  =  LINE  OUTPUT  DISK  UNIT  NUMBER. 

C  DATA  =  IMAGE  DATA 

C  I FLAG  =1  >>>  KL  TRANSFORM 

C  I FLAG  <>  1  >>>  ALL  OTHER  TRANSFORMS 

c******************************************************** 

C  USER  NOTE:  BE  CERTAIN  TO  DUPLICATE  ALL  COMMON  STATE- 
C  MENTS . 

c******************************************************** 
EXTERNAL  RTN 

REAL  DATA ( IB LKS Z  ,  N ) ,  A  ( 1 0  2  4 ) 

REAL *8  T  (16 , 16 ) 

INTEGER  OBSZ  , OTYPE 

COMMON/ CODPRM/  IR, IC, IBSZ , NTYPE, OBSZ , OTYPE, IBLKSZ , I FLAG 

COMMON/KLCOM/  T, IUNIT 

IF  (I FLAG  .EQ.  1)  REWIND  IUNIT 

DO  50  LINEl=l,M, IBLKSZ 

DO  20  1=1, IBLKSZ 

LINE  =  LINE1  +1-1 

CALL  DSKIO ( A, N, LINE, 1 , LIPT, NTYPE ) 

DO  10  J=1,N 
10  DATA ( I , J )  =  A { J) 

20  CONTINUE 

DO  30  IBLK=1,N, IBLKSZ 

£******************************************************** 

C  SOME  SPECIAL  I/O  HANDLING  OF  EIGENVECTOR  TRANS- 
C  FORMATION  IF  KL  TRANSFORMATION.  ( I FLAG  =  1) 

c******************************************************** 

IF ( (IR.EQ.-l.OR. IC.EQ.-l) .AND. (IFLAG.EQ.l) ) READ ( IUNIT) T 
CALL  RTN ( DATA ( 1, IBLK) , IBLKSZ , IBLKSZ, IR, IC, IBLKSZ , IBLKSZ ) 
30  IF( (IR.EQ.l.OR. IC.EQ.l) .AND. (IFLAG.EQ.l) ) WRITE ( IUNIT) T 
DO  50  1=1, IBLKSZ 
LINE  =  LINE1  +1-1 
DO  40  J=1 , N 
40  A ( J)  =  DATA ( I , J ) 

50  CALL  DSKIO (A, N, LINE, 0 , LOPT, OTYPE) 

RETURN 

END 
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nnoooooono 


*** 


Subroutine  to  Handle  Disk  Input  /  Output 


*** 


SUBROUTINE  DSKIO (A, N, LINE, RW , UNIT, TYPE) 
************************************************************** 

PURPOSE:  READS  AND  WRITES  INFORMATION  TO  AND  FROM  THE  DISK 

A  LINE  AT  A  TIME 

A=DATA  YOU  WANT  TOR  READ  OR  WRITE 
N=DATA  VECTOR  SIZE 
LINE=DATA  LINE  NUMBER 
RW = RE AD /WR I TE  FLAG 

UNIT=DISK  UNIT  NUMBER  FOR  I/O  OPERATION 
TYPE=DATA  TYPE  (E.  G.  INTEGER* 2,  REALM.  .  . ) 
************************************************************** 

INTEGER  TB(100)/100*0/ 

INTEGER  RW, UNIT, TYPE, C4 (1024) , TEST2/4HIN*2/ , TEST4/4HIN* 4/ 
INTEGER*2  B2(1024) 

REAL  A(N) 

IF (LINE. GT.TB( UNIT) )  GO  TO  10 
TB (UNIT)  =  0 
REWIND  UNIT 

10  TB (UNIT)  =  TB (UNIT)  +  1 
IF (RW. EQ . 0 )  GO  TO  500 
IF (TYPE . NE. TEST2 )  GO  TO  20 
READ  (UNIT)  (B2  ( I )  ,  1=1 ,  N) 

DO  15  1=1, N 
15  A ( I )  =  B2  ( I) 

RETURN 

20  IF(TYPE.NE.TEST4)  GO  TO  30 
READ (UNIT)  (C4 ( I ) , 1=1 , N) 

DO  25  1=1, N 
25  A ( I )  =  C4 ( I ) 

RETURN 

30  READ (UNIT)  (A(I),I=1,N) 

RETURN 

500  IF (TYPE. NE. TEST2)  GO  TO  520 
DO  515  1=1, N 
515  B2(I)  =  ( A ( I ) +.5) 

WRITE (UNIT)  (B2 ( I ) , 1=1 , N) 

RETURN 

520  IF (TYPE. NE. TEST4)  GO  TO  530 
DO  525  1=1, N 
525  C4 ( I)  =  A ( I ) 

WRITE (UNIT)  (C4 (I) , 1=1, N) 

RETURN 

530  WRITE ( UNIT)  (A(I),I=1,N) 

RETURN 

END 


*  ** 


Subroutine  for  Hadamard  Transform  *** 


SUBROUTINE  HADGEN (C, MC, NC, IR, IC,LR,LC) 

C 

C  COMPUTES  2  DIM.  HADAMARD  TRANSFORM 

C  OF  MATRIX. 

C  C  =  INPUT/OUTPUT  MATRIX 

C  MC=#OF  ROWS  OF  C 

C  NC=#OF  COLUMNS  OF  C 

C  IR  =  0  FOR  NO  ROW  TRANSFORM 

C  IC  =  0  FOR  NO  COLUMN  TRANSFORM 

C 

REAL  C (LR, LC) 

COMMON  /CORE/A (1024) 

C 

C  COLUMN  TRANSFORM 

C 

IF  (  IC  .EQ.  0  )  GO  TO  120 
CALL  POWER (  MC  ,  LMC  ) 

DO  110  N=1 , NC 
C 

C  FOR  EACH  COLUMN  TRANSFORM  ELEMENTS  1  THROUGH  MC 

C 

DO  111  IZAP=1 ,1024 
111  A ( IZAP) =0 . 

DO  112  M=1,MC 

A  (M)  =C  (M,  N) 

112  CONTINUE 

CALL  FHT  (LMC,  A) 

DO  115  M=1 , MC 

C  (M,  N)  =A  (M) 

115  CONTINUE 

110  CONTINUE 

C 

C  ROW  TRANSFORM 

C 

120  IF  (  IR  .EQ.  0  )  RETURN 

CALL  POWER (  NC  ,  LNC  ) 

DO  125  M=1 , MC 
C 

C  COPY  ROW  M  INTO  ARRAY 

C 

DO  130  N=1 , NC 

130  A ( N) =C (M, N) 

CALL  FHT (LNC, A) 

DO  140  N=1 , NC 

140  C (M, N) =A(N) 

125  CONTINUE 

RETURN 
END 
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***  Fast  Hadamard  Transform  *** 

SUBROUTINE  FHT  (M,X) 
DIMENSION  X (1) 

N=2**M 

NH=N/2 

NHP=3*NH 

NM=N-I 

DO  100  L=1 ,  M 
LTEST= (L/2) *2 
IF (LTEST. EQ . L)  GO  TO  200 
JJ=0 
JJN=NH 
JK2=N 
JK1=NM 
GO  TO  300 
200  JJ=N 

JJN=NHP 

JK2=0 

JK1=-1 

300  DO  100  J=1 , NH 

JJ=JJ+1 
JJN= JJN+1 
JKl=JKl+2 
JK2=JK2+2 

X(JKl)  -  X(JJ)  +  X(JJN) 

X  ( JK2)  =  X(JJ)  -  X(JJN) 
100  CONTINUE 

MTEST= (M/2) *2 
IF (MTEST. EQ. M)  GO  TO  400 
DO  500  1=1, N 
500  X (I) =X (N+I) 

400  RETURN 


***  Fast  Cosine  Transform  *** 


SUBROUTINE  FCT  ( INPUT,  ISIZE,  ITYPE) 

C  * 

C  THIS  SUBROUTINE  COMPUTES  THE  COSINE  TRANSFORM  OF  DATA  STORED 

C  IN  AN  INPUT  ARRAY  AND  RETURNS  THE  RESULTS  IN  THE  SAME  ARRAY. 

C  THERE  IS  A  CHECK  MADE  SO  THE  TABLE  IS  ONLY  CALCULATED  THE 
C  FIRST  TIME  A  CALL  IS  MADE,  OR  IF  ISIZE  CHANGES. 

C  MAXIMUM  SIZE  ARRAY  =  256  ELEMENTS 
C 

C  INPUT  =  INPUT/OUTPUT  ARRAY 

C  ISIZE  =  NUMBER  OF  ELEMENTS  IN  INPUT/OUTPUT  ARRAY 
C  ITYPE:  1=  FORWARD  TRANSFORM 

C  -1=  INVERSE  TRANSFORM 

C 

REAL  S  ( 12  8) 

REAL  COSINE (255) ,SINE(255) , INPUT(512) ,OUTPUT(1024) 

INTEGER  INV (128) ,MM(3) , OISIZE, OTYPE 
DATA  OTYPE/ 0/ 

DATA  MM/9,0,0/ 

DATA  ICHECK/0/ 

IF(ICHECK. NE.0. AND. ISIZE. EQ. OISIZE)  GOTO  100 

ICHECK=1 

OISIZE= ISIZE 

ISIZEF=ISIZE*4 

SQRT2=SQRT (2.0) 

CALL  POWER ( ISIZ E, MM (1 ) ) 

MM(1) =MM(1) +1 
C 

C  COSINE  TRANSFORM  TABLE  GENERATION 
C 

Y=3. 1415927/2. 0/FLOAT (ISIZE) 

CC=COS ( Y) 

SC=SIN ( Y) 

COSINE ( 1) =CC 
SINE ( 1 ) =SC 
JJ=ISIZE-2 
DO  950  1=1, JJ 

COSINE ( I +1 ) =COSINE ( I ) *CC  -  SINE(I)*SC 
950  SINE ( 1+1 ) =SINE ( I ) *CC  +  COSINE(I)*SC 
C  END  COSINE  TRANSFORM  TABLE  GENERATION 
C 

C  BRANCH  TO  DO  FORWARD  OR  INVERSE  TRANSFORM 
100  IF(ITYPE) 888,888,999 
C 

C  FORWARD  TRANSFORM 
C 

999  OUTPUT  (1)  =  INPUT  (1) 

OUTPUT ( 2)  =0.0 
DO  900  1=2, ISIZE 
OUTPUT  (1*2-1)  =INPUT(I) 

OUTPUT  ( ISIZ  EF+3- 1*2) =0.0 


non  non  non  non  ooonononooo 


***  Subroutine  for  Cosine  Transform  *** 

SUBROUTINE  COSINE (C, MC, NC, IR, IC,LR,LC) 

C= INPUT  /  OUTPUT  MATRIX 

MC=NUMBER  OF  ROWS  OF  C  TO  BE  TRANSFORMED 
NC=NUMBER  OF  COLUMNS  OF  C  TO  BE  TRANSFORMED 
IR=1,  FORWARD  ROW  TRANSFORM 
IR=0,  NO  ROW  TRANSFORM 
IR=~1 ,  INVERSE  ROW  TRANSFORM 
IC=1,  FORWARD  COLUMN  TRANSFORM 
IC=0,  NO  COLUMN  TRASFORM 
IC=-1,  INVERSE  COLUMN  TRANSFORM 

COMMON  /CORE/  A (1024) 

REAL  C (LR, LC) 

IF  (  7  .EQ.  0)  GO  TO  140 

REPEAT  FOR  EACH  COLUMN 

DO  110  N=1 , NC 

C ( 1 , N)  POINTS  TO  ROW  1  COLUMN  N 

DO  112  M=1 , MC 
A(M) =C(MrN) 

112  CONTINUE 

CALL  FCT ( A, MC, IC) 

DO  115  M=1,MC 
C (M,N) =A(M) 

115  CONTINUE 
110  CONTINUE 

140  IF  (IR  , EQ.  0)  RETURN 

FOR  EACH  ROW  COPY  TO  DUMMY  ARRAY 

DO  150  M=1 , MC 

COPY  ROW  M  INTO  A 

DO  130  N=1 , NC 
130  A(N) =C (M, N) 

CALL  FCT  (A,  NC,  IR) 

DO  145  N=1 , NC 
145  C (M, N) =A(N) 

150  CONTINUE 
RETURN 
END 
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***  Fast  Cosine  Transform  *** 

OUTPUT ( ISIZEF+4-I*2) =0.0 
900  OUTPUT (I *2) =0 .0 

OUTPUT (2*ISIZE+1) =0. 

OUTPUT (2*ISIZE+2) =0. 

CALL  HARM (OUTPUT, MM, INV,S,1, IFERR) 

IF ( IFERR. NE. 0 )  GOTO  998 
INPUT (l)=OUTPUT(l) *SQRT2 
DO  810  1=2  ,  ISIZE 

810  INPUT(I) =2.0* ( (OUTPUT (I *2-1 ) *COSINE(I-l) ) - 
*OUTPUT (1*2) *S INE ( 1-1 )  ) 

GOTO  1001 
C 

C  INVERSE  TRANSFORM 
C 

888  SUM=  INPUT  (1)/SQRT2 

CT=SUM* (1.0-1. 0/SQRT2 ) /FLOAT (ISIZE) 

OUTPUT  (1)  =  INPUT  (1) 

OUTPUT  (2)  =0.0 
DO  910  1=2  r ISIZ  E 

OUTPUT ( 1*2-1 ) =INPUT ( I ) *COSINE(I-l) 
OUTPUT(ISIZEF+3-I*2) = INPUT ( I ) *COSINE ( 1-1 ) 
OUTPUT ( ISIZEF+4-I*2) =  INPUT ( I ) *SINE ( 1-1 ) 

910  OUTPUT (I*2)=-(INPUT(I)  *SINE(I-1)  ) 

OUTPUT ( 2*IS IZE+1 ) =0. 

OUTPUT (2*ISIZE+2) =0 . 

CALL  HARM (OUTPUT, MM, INV, S, -1 , IFERR) 

C  IBM  FFT  ROUTINE  (3D.) 

IF ( IFERR. EQ.0)  GOTO  960 
998  IF (OTYPE. NE. 0 ) WRITE (6,612)  IFERR 

612  FORMAT ('  ' , ' HARM  ERROR' , 1  *  ERROR=  ’,12) 

STOP  777 
960  CONTINUE 

DO  811  1=1, ISIZE 

811  INPUT ( I) =OUTPUT( 1*2-1 ) +CT 
1001  CONTINUE 

RETURN 

END 


***  Subroutine  to  Determine  Power  *** 


SUBROUTINE  POWER(JJ, IPOWER) 

C**  THIS  SUBROUTINE  INPUTS  A  VALUE  ' J J '  AND  RETURNS 
C**  ITS  POWER  OF  2  -  'IPOWER'. 

C**  WHERE: 

C**  IPOWER 

C**  JJ  =  2 

C** 

INTEGER  OTYPE 
DATA  OTYPE/ 0/ 

DATA  IOUT/6/ 

C  IOUT  IS  THE  PRINTER  UNIT  NUMBER. 

NOJJ 
IPOWER=0 
10  I=NC/2 

NC=I 

IPOWER=IPOWER+l 
IF(I.GT.l)  GOTO  10 
IF(I.EQ.l)  GOTO  20 
IF (OTYPE. NE. 0) WRITE (IOUT, 1) 

1  FORMAT ( '  ','****  POWER  ERROR  ****  INPUT  IS  NOT  A  POWER  OF  2 

STOP  777 
20  RETURN 


I 

: 

•« 

I 

j 

i 


nnnonnoonooonoonnnnoono 


Main  Driver  Program  for  all  ASPC  Simulations  *** 

INTEGER  IN2/ 4HIN*2/ , IN4/ 4HIN*4/ , RL/4HREAL/ 

REAL *4  XIN(256,256) , COEFF ( 7 , 256 ) 

REAL *4  QTH ( 50) ,QL(50) 

REAL  A (256) 


AH PC/AS PC  MAIN  DRIVER  PROGRAM 

EXPLANATIONS  OF  VARIABLES: 

XIN  -  MAIN  I/O  IMAGE  MATRIX 

COEFF  -  COEFFICIENT  STORAGE  MATRIX 

JBLKSZ  -  TRANSFORMATION  BLOCKSIZE 

ISZ  -  IMAGE  ROW  SIZE 

JSZ  -  IMAGE  COLUMN  SIZE 

NUMSMP  -  NUMBER  SAMPLES/LINE  (SAME  AS  JSZ) 

NUMFR  -  NUMBER  SAMPLE  LINES/FRAME  (SAME  AS  ISZ ] 
QTH  -  QUANTIZER  THRESHOLDS 

QL  -  QUANTIZER  LEVELS 

SUBROUTINES: 

DSKIO  -  DISK  INPUT  /  OUTPUT 
MAXQTZ  -  MAX  QUANTIZER  DESIGN 
COMPMC  -  ADPCM  DATA  COMPRESSION 


NUMFR=256 

NUMSMP=256 

NCOEFF=3 

ISZ=256 

JSZ=256 

JBLKSZ=16 

COMMON  /A/  JBLKSZ 

DO  10  1=1, ISZ 

CALL  DSKIO (A, JSZ, 1,1,1, RL) 

DO  10  J=1,JSZ 

10  XIN ( I , J) =A ( J) 

DO  11  1=1,12 

11  WRITE (6,1001) (XIN ( I , J) , J=1 , 16 ) 

1001  FORMAT ('  ' ,16 (F6.0,1X) ) 

NBIT=1 

CALL  MAXQTZ (QTH, QL,NLEVEL,NB IT) 

CALL  COMPMC (XIN, NUMFR, NUMSMP, NCOEFF, COEFF, QTH , QL, NLEVEL) 
WRITE(6 ,20) 

20  FORMAT ('  ******  MADE  IT  BACK  THROUGH  SUBS  *****★') 

DO  30  1=1,12 

30  WRITE (6,1001) (XIN ( I , J) , J=1 , 16 ) 

DO  40  1=1, ISZ 
DO  35  J=1,JSZ 
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***  Main  Driver  Program  for  all  A S?C  Simulations  *** 


35  A ( J) =XIN ( I ,  J) 

Aa  <£LL  D^KIO  (A,  JSZ  / 1, 0 , 2 ,  RL) 

40  CONTINUE 
WRITE (6,2) 

2  FORMAT ('0***  JOB  FINISHED  ***') 
STOP  ' 

END 
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non  n n n n n n n o o n n 


***  Subroutine  for  ADPCM  using  ARMA  Model 


*  *  * 


SUBROUTINE  COMPMC (XIN, NF, NSF, N, CM, QTH , QL, NLEVEL) 

REAL  A(10) , XIN (NF, NSF) ,  R(  7) , XV ( 256 ) , XTR ( 256 , 256 ) 

REAL  CM ( 7 , 2 5 6 ) ,QL(1) ,QTH(1) ,E2V(256) ,VAREST(256) ,XMN(256) 

REAL  ACV (256 ) ,AC(10) ,PACV(10) , WORK (10) 

REAL  P ( 10 ) fWW(150) ,RESMN(256) 

REAL  B(10,10) , PAR (10) 

REAL  QML (1024) ,QMTH (1024) 

INTEGER  NARPS ( 256 ) 

EXPLANATION  OF  CALL  VARIABLES: 

XIN  -  THE  INPUT  DATA  MATRIX  AND  RESIDUAL  OUTPUT  MATRIX 

NF  -  NUMBER  OF  FRAMES  (NUMBER  OF  LINES  IN  THE  INPUT  MATRIX) 

NSF  -  NUMBER  OF  SAMPLES  PER  FRAME  (SAMPLES  PER  LINE) 

N  -  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED 

CM  -  MATRIX  THAT  CONTAINS  THE  PREDICTOR  COEFFICIENT  VECTORS 


SET  UP  THE  CONSTANTS  AND  INITIAL  VALUES 

CALL  MAXQTZ (QMTH, QML, MLEVEL, 6 ) 

XNSF  =  NSF 
COMMON  /A/  JBLKSZ 
XAVBLK=JBLKSZ  *12  8. 

IP=N 

IQ=0 

NHOLD=N 

IN  THE  DO  160  LOOP,  I  IS  THE  LINE  NUMBER  (1  -  NF) 

VSME=0 . 

VSMSE=0. 

RSME=0  . 

RSMSE=0 . 

NLEV=2*NLEVEL 
DO  160  1=1, NF 
N=NHOLD 
SUM=0 . 

SUMSX=0. 

DO  501  J=1 , NSF 
501  XV(J) =XIN(I, J) 

DO  510  J=1 , NSF 
SU  M=  SU  M+XV  ( J ) 

SUMSX=SUMSX+XV  ( J)  *XV  ( J) 

510  CONTINUE 

XMN(I) =SUM/NSF 

IF(MOD(I-l, JBLKSZ)  .NE.  0)GO  TO  373 
XMN ( I ) =2750 . 

GO  TO  374 

373  XMN ( I ) =0 . 

374  CONTINUE 


o  o 


***  Subroutine  for  ADPCM  using  ARMA  Model 


*  *  * 


321  CALL  FTAUTO {XV , NSF, N, N, 7 , XMN ( I )  , ACV ( 1 ) , ACV ( 2 )  , AC, PA  CV,  WORK) 

IP=N 

IQ=0 

CALL  FT ARPS (ACV / XBAR, IP, IQ, A, PMAC,WW, IER) 

IF (N  .EQ.  2) GO  TO  432 
IF  (IER  .NE.  129)  GO  TO  432 
N=  N-l 

WRITE ( 6 , 636) N 

636  FORMAT ('  - >  N  MADE  TO:  ',13) 

GO  TO  321 

432  CALL  PARSBL(A, N, PAR, B) 

DO  335  KK=1,N 

CM (KK, I) -UNI FRM (-1.75,1. 75, MLEV EL, A(KK) ) 

335  CONTINUE 

DO  700  *11=1, N 
700  R ( I I ) =XMN ( I ) 

DO  710  11=1, NSF 
PRED=0 . 

DO  705  J=1,N 

705  PRED=CM ( J , I ) *R { J ) +PRED 
PRED=PRED+PMAC 
E2V(II) =XV (II) -PRED 

DO  706  K=2,N 

706  R(N+2-K) =R (N+l-K) 

R ( 1) =XV (II) 

710  CONTINUE 
SUME=0 . 

SUMSQE=0 . 

DO  155  M=1 , NSF 
SUME=SUME+E2 V(M) 

155  SUMSQE=SUMSQE+ ( E2 V(M) *E2V(M) ) 

VAREST ( I ) = ( SUMSQE- ( (SUME*SUME) /NSF) ) / (NSF-1 . ) 

VSME=VSME+VAREST ( I ) 

VSMSE=VSMSE+VAREST ( I ) *VAREST ( I ) 

RESMN ( I ) =SUME/NSF 

WRITE (6,919)1, RESMN ( I ) ,VAREST(I) , (A(K1) ,K1=1,N) , (CM(K2,I) ,K2=1,N) 
919  FORMAT ( 1  ' , 14 , 5 ( IX, E10 . 3 ) , / , 3 ( IX, E10 . 3 ) ) 

VAREST ( I ) =UNIFRM ( 0 . , 40000 512 , VAREST ( I ) ) 

IF (MOD (1-1 , JBLKSZ )  .EQ.  0 ) RESMN ( I ) =RESMN ( I ) /JBLKSZ 
RSME=RSMN+RESMN ( I ) 

RSMSE=RSMSE+RESMN { I ) *RESMN ( I ) 

RESMN { I ) =QNTZ ( RESMN ( I) , MLEV EL , QML , QMTH ,150.) 

IF(MOD(I-l, JBLKSZ)  .EQ.  0) RESMN ( I ) =RESMN ( I ) *JBLKSZ 
DO  200  J=1 , N 
200  R ( J )  =  XMN ( I ) 

DO  250  J=1 , NSF 
S=XV(J) 

RT2=0 .0 

NOW  GET  XMITTER  RESIDUAL 
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***  Subroutine  for  ADPCM  using  ARMA  Model  *** 

C 

DO  230  K=1 , N 

230  RT2=RT2  +  CM{K,I}*R(K) 

E2  =  S  -  RT2  -  RESMN(I) 

EQ2=QNTZ (E2 , NLEVEL, QL, QTH, VAREST ( I ) ) 

XTR ( I ,  J) =EQ2 
C 
C 

C  SHIFT  R-REGISTER 

C 

C 

DO  240  K=2,N 
240  R (N+2-K)  =  R(N+1-K) 

R ( 1)  =  RT2  +EQ2  +RESMN ( I ) 

250  CONTINUE 
NARPS ( I ) =N 
160  CONTINUE 

VMN=VSME/256. 

WAR=  ( VSMSE-  (  (VSME*VSME)/256.)  )/255. 

RMN=RSME/256. 

RVAR= (RSMSE- ( (RSME*RSME) /256 . ) )/255. 

WRITE  ( 6 , 920 )  VMN ,  WAR,  RMN ,  RVAR 

920  FORMAT ( 1  VMN= • , El 5 . 5 , '  W  1  ,E15.5,'  RMN= ' , El 5 . 5 , 1  RVAR=',E15.5 

c 

C  RECEIVER  SIMULATION  SECTION 
C 

- - 

c 

C  SET  R  -  REGISTER  TO  MEAN  VALUE 
C  FOR  THE  CURRENT  LINE. 

C 

C 

DO  210  I  =  1 , NF 
N= NARPS ( I ) 

DO  170  J  =  1 , N 
170  R( J) =XMN ( I ) 

DO  400  J  =  1 , NSF 
RV2  =0.0 
DO  180  K  =  1,N 
180  RV2  =  RV2  +  CM  (K,  I)  *R  (K) 

RN  =  XTR ( I , J)  +  RV2  +  RESMN(I) 

XIN ( I , J) =RN 

C - 

C 

C  SHIFT  R  -  REGISTER. 

C 

C 

DO  190  K  =  2 ,  N 
190  R (N+2-K)  =  R(N+1-K) 
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***  Subroutine  for  ADPCM  using  Kalman  Filter 


*** 


SUBROUTINE  COMPMC (XIN, NF, NSF , N, CM, QTH , QL, NLEVEL) 

REAL  A (10) ,V1 (10) , XIN (256, 2 56) ,R(10) , XV (256) ,G(10) ,V(10,10) 

REAL  XTR(256 ,256) 

REAL  CM ( 7 , 2 5 6 ) , QL ( 1 ) , QTH ( 1 ) ,E2V(256)  ,VAREST(256)  ,XMEAN(256) 

REAL  RESMN ( 256 ) 

REAL  QML ( 1024) ,QMTH( 1024) 

REAL  B(10 ,10) , PAR (10) 

EXPLANATION  OF  CALL  VARIABLES: 

XIN  -  THE  INPUT  DATA  MATRIX  AND  RESIDUAL  OUTPUT  MATRIX 

NF  -  NUMBER  OF  FRAMES  (NUMBER  OF  LINES  IN  THE  INPUT  MATRIX) 

NSF  -  NUMBER  OF  SAMPLES  PER  FRAME  (SAMPLES  PER  LINE) 

N  -  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED 

CM  -  MATRIX  THAT  CONTAINS  THE  PREDICTOR  COEFFICIENT  VECTORS 

DEFINITION  OF  VARIABLE  TERMS: 

V  -  THE  ERROR  COVARIANCE  MATRIX 

A  -  THE  PREDICTOR  COEFFICIENT  VECTOR 

VAR I  -  INITIAL  VALUE  FOR  THE  ERROR  COVARIANCE  MATRIX 

W  -  VARIANCE  OFFSET 

XV  -  THE  INPUT  LINE  TEMPORARY  VECTOR 

G  -  THE  GAIN  VECTOR 

R  -  THE  PAST  VALUE  VECTOR 

E  -  THE  ERROR  OR  RESIDUAL  TERM 

SET  UP  THE  CONSTANTS  AND  INITIAL  VALUES 

DATA  V/ 100*0.0/, A/ 1.0, -.5, -.2, .3, .4, -.5/ 

CALL  MAXQTZ ( QMTH , QML , MLEVEL, 6 ) 

VLO=  9999999. 

VHIGH=-1. 

XNSF  *  NSF 
COMMON  /A/  JBLKSZ 
XAVBLK* JBLKSZ *128. 

NSFM1=NSF-1 
W  =  1.0 

IN  THE  DO  160  LOOP,  I  IS  THE  LINE  NUMBER  (1  -  NF) 

DO  160  1=1, NF 
VAR 1=1 00 . 

DO  10  J=1 ,  N 
DO  10  K=1 ,  N 
V(J,K)  *  0.0 

IF(J.EQ.K)  V(J,K)  =  VARI 
10  CONTINUE 

SET  UP  THE  INPUT  VECTOR  AND  THE  PAST  VALUE  VECTOR 


% 


SUMX=0 
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***  Subroutine  for  ADPCM  using  Kalman  Filter  *** 

DO  20  J=1 , NSF 
SUMX=SUMX+XIN(I, J) 

20  XV(J)  =  X1N  ( I ,  J) 

XMEAN ( I ) =SUMX/XNSF 
WRITE (6,2000) I, XMEAN( I) 

2000  FORMAT ( 1  LINE:  ' , 15 ,2X, E15 .5) 

IF (MOD (1-1 , JBLKSZ )  .NE.  0)GO  TO  373 
XMEAN (I) =2750. 

GO  TO  374 

373  XMEAN ( I) =0  . 

374  DO  30  J=1 , N 
30  R ( J)  =  XMEAN ( I ) 

IDENTIFICATION  LOOP  (IDENTIFY  THE  PREDICTOR  COEFFICIENTS) 


DO  100  J=1 , NSF 
S  =  XV(J) 

S2  =  0.0 
RT  =  0.0 
DO  40  K=1,N 
VI (K)  =  0. 

DO  40  L=1 ,  N 

40  VI (K)  =  Vl (K)  +  V(K, L) *R (L) 

DO  50  K=1 , N 

50  S2  *  S2  +  R  (K)  *Vl  (K) 

DO  60  K=1 , N 

G  (K)  =  Vl  (K)  /  ( VARI  +  S2) 

60  RT  =  RT  +  A  (K)  *R  (K) 

DO  70  K=1 , N 
DO  70  L=1 , K 

V(K,L)  =  V(K, L)  -  G ( K) *V1 ( L) 

70  V(L,K)  =  V(K, L) 

E  =  S  -  RT 
DO  80  K=1 , N 

80  A (K)  =  A (K)  +  G (K)  *E 

SHIFT  THE  PAST  VALUE  VECTOR 
DO  90  L=2 , N 

90  R  (N+2-L)  =  R(N-*-l-L) 

R ( 1)  =  RT+E 

IF ( I  .EQ.  225) WRITE (6 ,1001) (G(K1) ,K1=1,N) 
1001  FORMAT ( '  G= ' , 3 ( IX , E15 .5 ) ) 

100  CONTINUE 

CALL  PARSBL(A, N, PAR, B) 

DO  335  KK=1 , N 

CM (KK, I) =UNIFRM ( -1.75,1. 75, 256, A(KK) ) 

335  CONTINUE 


***  Subroutine  for  ADPCM  using  Kalman  Filter  *** 

C  RESET  R-REGISTER  TO  INITIAL  SIGNAL  VALUE 

C 

DO  120  J=1  ,N 
120  R ( J)  =  XMEAN(I) 

DO  150  J=1 , NSF 
S  =  XV (J) 

RT2  =0.0 
C 

C  DETERMINE  RESIDUAL  SIGNAL 

C 

C  RT2  IS  AT  THE  TRANSMITTER  STAGE  WHERE: 

C  RT2  =  RT2  +  AO  (K)  *RV  (K) 

C 

DO  130  K=1 , N 

130  RT2  =  RT2  +  CM(K,I)*R(K) 

E2V( J) =S-RT2 
C 

C  SHIFT  R-REGISTER 

C 

DO  140  K=2 , N 
140  R(N+2-K)  =  R (N+l-K) 

R  ( 1)  =  RT2+E2V{  J) 

150  CONTINUE 
SUME=0 . 

SUMSQE=0 . 

DO  155  M=1 , NSF 
SUME=SUME+E2V(M) 

155  SUMSQE=SUMSQE+ (E2V(M) *E2V(M) ) 

VAREST ( I ) = ( SUMSQE- ( SUME*SUME/NSF) ) /NSFM1 
IF ( VAREST ( I )  .LT.  VLO) VLO=VAREST ( I ) 

IF (VAREST ( I )  .GT.  VHIGH) VHIGH= VAREST ( I ) 

VSME=VSME+VAREST ( I ) 

VSMSE=VSMSE+ VAREST ( I ) *VAREST ( I ) 

RESMN ( I ) =SUME/XNSF 

IF (MOD ( I — 1 , JBLKSZ )  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) /JBLKSZ 

WRITE (  6 , 91  9)  I , RESMN ( I )  , VAREST (I)  , (A(K1)  ,K1=1,N)  , (CM(K2,I)  ,K2=1,N) 
919  FORMAT { '  ' ,14,5 (IX, E10 .3) ,/,3 (IX, E10.3) ) 

VAREST(I)=UNIFRM(0. ,40000. ,512 , VAREST (I) ) 

RSME= RSME+RESMN ( I ) 

RSMSE=RSMSE+RESMN ( I ) *RESMN ( I ) 

RESMN ( I ) =QNTZ (RESMN ( I ) , MLEVEL , QML , QMTH ,150.) 

IF (MOD (1-1, JBLKSZ)  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) *JBLKSZ 
DO  200  J=1 ,N 
200  R{ J)  =  XMEAN(I) 

DO  250  J=1 , NSF 
S=XV(J) 

RT2=0 .0 
C 

C  NOW  GET  XMITTER  RESIDUAL 
C 
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***  Subroutine  for  AD PCM  using  Kalman  Filter  *** 

DO  230  K=1 , N 

230  RT2=RT2  +  CM(K,I)*R(K) 

E2  =  S  -  RT2  -  RESMN(I) 

EQ2=QNTZ (E2,NLEVEL,QL,QTH, VAREST(I) ) 

XTR ( I , J) =EQ2 
C 
C 

C  SHIFT  R- REGISTER 

C 

C 

DO  240  K=2 , N 
240  R ( N+2-K)  =  R{N+1-K) 

R(l)  =  RT2+EQ2+RESMN (I) 

250  CONTINUE 
160  CONTINUE 

VMN=VSME/256 . 

WAR=  ( VSMSE-  (  (VSME*VSME)  /256  .  )  )  /255  . 

RMN=RS ME/256. 

RVAR= (RSMSE- ( (RSME*RSME) /256 . ) ) /255. 

WRITE  (6,920)  VMN ,  WAR,  RMN ,  RVAR 

920  FORMAT  (  '  VMN=  '  ,  El  2 . 5  ,  '  WAR=  •  ,  El  2 . 5  ,  '  RMN=  1  ,  El  2 . 5  ,  '  RVAR=  1 
WRITE ( 6 , 921 ) VLO, VHIGH 

921  FORMAT ('  VLO= ' , El  5 . 5 , '  VHIGH= ' , El 5 . 5 ) 

- - 

C 

C  RECEIVER  SIMULATION  SECTION 
C 

- - 

C 

C  SET  R  -  REGISTER  TO  MEAN  VALUE 
C  FOR  THE  CURRENT  LINE. 

C 

C 

DO  210  I  =  1 , NF 
DO  170  J  =  1,N 
170  R ( J) =XMEAN ( I ) 

DO  400  J  =  1 ,NSF 
S  =  XIN ( I , J) 

RV2  =0.0 
DO  180  K  =  1,N 
180  RV2  =  RV2  +  CM  (K,  I)  *R  (K) 

RN  =  XTR ( I , J)  +  RV2  +  RESMN(I) 

XIN ( I , J)  =  RN 

- - 

C 

C  SHIFT  R  -  REGISTER. 
r 

C 

DO  190  K  =  2 ,  N 
190  R (N+2-K)  =  R (N+l-K) 
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***  Subroutine  for  AD PCM  using  Stoch.  Approx. 


*  ★  ★ 


SUBROUTINE  COMPMC (XIN, NF, NSF, N, CM, QTH , QL, NLEVEL) 

REAL  A(10) , XIN ( NF , NSF) ,R(7) ,XV(256) , XTR ( 256 , 2 56 ) 

REAL  CM ( 7 , NF) ,QL(1) ,QTH(1) ,E2V(256) ,VAREST(256) ,XMN(256) 

REAL  RESMN ( 256 ) 

REAL  B ( 10 , 10 ) , PAR (10) 

REAL  QML (1024) , QMTH (1024) 

EXPLANATION  OF  CALL  VARIABLES: 

XIN  -  THE  INPUT  DATA  MATRIX  AND  RESIDUAL  OUTPUT  MATRIX 

NF  -  NUMBER  OF  FRAMES  (NUMBER  OF  LINES  IN  THE  INPUT  MATRIX) 

NSF  -  NUMBER  OF  SAMPLES  PER  FRAME  (SAMPLES  PER  LINE) 

N  -  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED 

CM  -  MATRIX  THAT  CONTAINS  THE  PREDICTOR  COEFFICIENT  VECTORS 


IN  THE  DO  160  LOOP,  I  IS  THE  LINE  NUMBER  (1  -  NF) 

CALL  MAXQTZ (QMTH , QML , MLEVEL , 6 ) 

COMMON  /A/  JBLKSZ 
XAVBLK=JBLKSZ  *256 . 

DO  160  1=1, NF 
DO  501  J=1 , NSF 
501  XV ( J) =XIN ( I , J) 

CALL  STAPRX (XV , NSF, N, A, XMN ( I ) ) 

IF (MOD (1-1, JBLKSZ)  .NE.  0)GO  TO  373 
XMN ( I ) =2750  . 

GO  TO  374 

373  XMN ( I ) =0  . 

374  CALL  PARSBL ( A, N, PAR, B) 

DO  335  KK=1 , N 

CM(KK,I)=UNIFRM(-1.75,1.75,256,A(KK) ) 

335  CONTINUE 

DO  700  11=1, N 
700  R ( I I ) =XMN ( I ) 

DO  710  11=1, NSF 
PRED=0 . 

DO  705  J=1 , N 

705  PRED=CM(J, I) *R(J) +PRED 
E2 V( II ) =XV (II) -PRED 

DO  706  K=2 , N 

706  R(N+2-K) =R(N+1-K) 

R ( 1 ) =XV (II) 

710  CONTINUE 
SUME=0 . 

SUMSQE=0 . 

DO  155  M=1 , NSF 
SUME=SUME+E2V(M) 

155  SUMSQE=SUMSQE+(E2V(M) *E2V(M) ) 

VAREST ( I ) = ( SUMSQE- ( (SUME*SUME) /NSF) )/(NSF-l.) 


nonnooonon  r.nonn  non 


***  Subroutine  for  ADPCM  using  Stoch.  Approx.  *** 

RESMN { I ) =SUME/NSF 

IF (MOD ( I — 1 , JBLKSZ )  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) /JBLKSZ 

WRITE (6, 91 9) I, RESMN (I) ,VAREST(I) , (A(K1) ,K1=1,N) , (CM (K2 , I 
919  FORMAT ( '  ' , 14 , 5 ( IX, E10 .3 ) , /, 3 (IX, E10 .3 ) ) 
VAREST(I)=UNIFRM(0. ,100000. ,512, VAREST ( I) ) 

RESMN ( I ) =QNTZ ( RESMN ( I ) , MLEVEL ,  QML , QMTH ,150.) 

IF (MOD (1-1, JBLKSZ)  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) * JBLKSZ 
DO  200  J=1 , N 
200  R ( J)  =  XMN(I) 

DO  250  J=1 , NSF 
S=XV(J) 

RT2=0  .0 

NOW  GET  XMITTER  RESIDUAL 
DO  230  K=1 , N 

230  RT2=RT2  +  CM(K,I)*R(K) 

E2  =  S  -  RT2  -  RESMN (I) 

EQ2=QNTZ ( E2 , NLEVEL , QL, QTH , VAREST (I) ) 

XTR ( I , J) =EQ2 


SHIFT  R- REGISTER 


DO  240  K=2 , N 
240  R (N+2-K)  =  R(N+1-K) 

R ( 1)  =  RT2  +EQ2  +RESMN ( I ) 
250  CONTINUE 
160  CONTINUE 


RECEIVER  SIMULATION  SECTION 


SET  R  -  REGISTER  TO  MEAN  VALUE 
FOR  THE  CURRENT  LINE. 


DO  210  I  =  1 , NF 
DO  170  J  =  1 , N 
170  R ( J) =XMN ( I ) 

DO  400  J  *  1 , NSF 
RV2  =  0.0 
DO  180  K  =  1 ,  N 
180  RV2  =  RV2  +  CM  (K,  I)  *R  (K) 

RN  =  XTR ( I , J)  +  RV2  +  RESMN (I) 
XIN ( I , J) =RN 


)  , K2=l , N) 


C 


***  Subroutine  for  ADPCM  using  Stoch.  Approx 


C  SHIFT  R  -  REGISTER. 

C 

C 

DO  190  K  =  2 ,  N 
190  R (N+2-K)  =  R(N+1-K) 
R  ( 1)  =  RN 
400  CONTINUE 
210  CONTINUE 
RETURN 


***  Subroutine  to  Provide  PARCOR  Stabilization  *** 

SUBROUTINE  PARSBL { A, N, PAR, B) 

- - 

C 

C  COEFFICIENT  PARCOR  STABILIZATION 
C  SUBROUTINE 
C 
C 

REAL  A{1) ,  B  ( 1 0 , 1 0 ) , R ( 1 3 ) , PAR ( 1 ) 

10  DO  100  1=1 , N 
100  B ( I , N) =-A ( I ) 

IF (ABS ( A (N) )  .GT.  0.97)GO  TO  500 
PAR (N) =-A ( N) 

NM1=N-1 

DO  300  1=1, NM1 

K=N+1-I 

KM=K-1 

G=1.0/(1.0-B(K,K) *B(K,K) ) 

DO  200  KK=1 , KM 

200  B(KK, KM) = (B(KK, K) -B (K , K) *B (K-KK, K) ) *G 
IF (ABS ( B ( KM, KM) )  .GT.  0.97)GOTO500 
PAR (KM) =B (KM, KM) 

300  CONTINUE 
RETURN 
500  GG=1 . 0 

DO  510  1=1, N 
GG=0 . 97*GG 
510  A(I) =GG*A(I) 

GO  TO  10 
END 
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***  Subroutine  for  ADPCM  /  Kalman  filter  /  Zeroing  *** 

SUBROUTINE  COMPMC (XIN,  NF,  NSF,  N,  CM, QTH , QL, NLEVEL) 

REAL  A(10) ,Vl(10) , XIN (256,256) ,R(10) , XV (256) ,G(10) , V( 10 , 1 0 ) 

REAL  XTR (256 , 256 ) 

REAL  CM ( 7 ,256 ) ,QL(1) ,QTH(1) ,E2V(256) ,VAREST(256) ,XMEAN(256) 

REAL  RESMN ( 256 ) 

REAL  QML (1024) , QMTH (1024) 

REAL  B( 10 ,10) , PAR ( 10 ) 

EXPLANATION  OF  CALL  VARIABLES: 

XIN  -  THE  INPUT  DATA  MATRIX  AND  RESIDUAL  OUTPUT  MATRIX 

NF  -  NUMBER  OF  FRAMES  (NUMBER  OF  LINES  IN  THE  INPUT  MATRIX) 

NSF  -  NUMBER  OF  SAMPLES  PER  FRAME  (SAMPLES  PER  LINE) 

N  -  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED 

CM  -  MATRIX  THAT  CONTAINS  THE  PREDICTOR  COEFFICIENT  VECTORS 

DEFINITION  OF  VARIABLE  TERMS: 

V  -  THE  ERROR  COVARIANCE  MATRIX 

A  -  THE  PREDICTOR  COEFFICIENT  VECTOR 

VAR I  ~  INITIAL  VALUE  FOR  THE  ERROR  COVARIANCE  MATRIX 
W  -  VARIANCE  OFFSET 

XV  -  THE  INPUT  LINE  TEMPORARY  VECTOR 

G  -  THE  GAIN  VECTOR 

R  -  THE  PAST  VALUE  VECTOR 

E  -  THE  ERROR  OR  RESIDUAL  TERM 


SET  UP  THE  CONSTANTS  AND  INITIAL  VALUES 

DATA  V/ 100*0.0/ , A/1 .0,-.5,-.2,.3, .4, -.5/ 

LIMCNT=2 
IB EAT=1 

CALL  MAXQTZ (QMTH , QML, MLEVEL, 6 ) 

VLO=  9999999. 

VHIGH=-1 . 

XNSF  =  NSF 
COMMON  /A/  JBLKSZ 
XAVBLK= JBLKS  Z *12  8. 

NSFM1=NSF-1 
W  =  1.0 
VARI  =  100.0 

IN  THE  DO  160  LOOP,  I  IS  THE  LINE  NUMBER  (1  -  NF) 

DO  160  1=1, NF 
DO  10  J=1 , N 
DO  10  K=1 , N 
V ( J, K)  =  0.0 

IF(J.EQ.K)  V(J,K)  =  VARI 
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***  Subroutine  for  ADPCM  /  Kalman  filter  /  Zeroing  *** 

10  CONTINUE 

SET  UP  THE  INPUT  VECTOR  AND  THE  PAST  VALUE  VECTOR 

SUMX=0 . 

DO  20  J=1 , NSF 
SUMX=SUMX+XIN(I, J) 

20  XV(J)  =  XIN{I,J) 

IF (MOD ( I — 1 , JBLKSZ )  .NE.  0)GO  TO  373 
XMEAN(I) =2750. 

GO  TO  374 

373  XMEAN { I ) =0  . 

374  DO  30  J=1 , N 
30  R ( J)  =  XMEAN (I) 

IDENTIFICATION  LOOP  (IDENTIFY  THE  PREDICTOR  COEFFICIENTS) 


DO  100  J=1,NSF 
S  =  XV  (J) 

S2  =  0.0 
RT  =  0.0 
DO  40  K=1,N 
VI (K)  =  0. 

DO  40  L=1,N 

40  Vl (K)  =  Vl (K)  +  V(K,L)*R(L) 

DO  50  K=1 , N 

50  S2  =  S2  +  R (K) *Vl ( K) 

DO  60  K=1 , N 

G  (K)  =  Vl  (K )  /  (VAR I  +  S2) 

60  RT  =  RT  +  A (K ) *R (K) 

DO  70  K=1,N 
DO  70  L=1 , K 

V(K,L)  =  V(Kf  L)  -  G (K) *Vl ( L) 

70  V(L,K)  =  V(K,L) 

E  =  S  -  RT 
DO  80  K=1 ,  N 

80  A (K)  =  A (K)  +  G(K)*E 

SHIFT  THE  PAST  VALUE  VECTOR 
DO  90  L=2 , N 

90  R (N+2-L)  =  R(N+1-L) 

R ( 1 )  =  RT+E 
100  CONTINUE 

CALL  PARSB L ( A , N f  PAR , B ) 

DO  335  KK=1 , N 

CM(KK,  I)=UNIFRM(-1.75,1.75,'256,A(KK)  ) 
335  CONTINUE 


non  non  o  n  o  n  n  o  no 


***  Subroutine  for  ADPCM  /  Kalman  filter  /  Zeroing  *** 

RESET  R-REGISTER  TO  INITIAL  SIGNAL  VALUE 

DO  120  J=1 , N 
120  R(J)  =  XMEAN(I) 

DO  150  J=1 , NSF 
S  =  XV (J) 

RT2  =0.0 

DETERMINE  RESIDUAL  SIGNAL 

RT2  IS  AT  THE  TRANSMITTER  STAGE  WHERE: 

RT2  =  RT2  +  AO  (K)  *RV  (K) 

DO  130  K=1 , N 

130  RT2  =  RT2  +  CM(K,I)*R(K) 

E2V( J) =S-RT2 

SHIFT  R-REGISTER 

DO  140  K=2,N 
140  R (N+2-K)  =  R(N+1— K) 

R  ( 1)  =  RT2+E2V  { J) 

150  CONTINUE 
SUME=0 . 

SUMSQE=0 . 

DO  155  M=1,NSF 
SUME=SUME+E2V (M) 

155  SUMSQE=SUMSQE+(E2V(M)*E2V(M) ) 

VAREST(I) = (SUMSQE- (SUME*SUME/NSF) )/NSFMl 
IF ( VAREST ( I )  .LT.  VLO) VLO=VAREST ( I ) 

IF (VAREST(I)  ,GT.  VHIGH) VH IGH= VAREST ( I ) 

VSME=VSME+VAREST ( I ) 

VSMSE=VSMSE+VAREST ( I ) *VAREST(I) 

RESMN ( I ) =SUME/NSF 

IF (MOD ( 1-1 , JBLKSZ )  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) /JBLKSZ 

WRITE (6, 91 9) I, RESMN { I) f VAREST ( I ) , (A (Kl) ,K1=1,N) , (CM (K2, I) ,K2 
919  FORMAT ( '  ' , 14 ,5 (IX, E10 .3) ,/,3 (IX, E10 .3) ) 

VAREST (I) =UNIFRM (0. ,40000 . ,512, VAREST (I) ) 

RSME=RSME+RESMN ( I ) 

RSMS  E=  RSMS  E+ RESMN ( I ) *RESMN ( I ) 

RESMN ( I ) =QNTZ (RESMN (I) , MLEVEL, QML, QMTH ,150.) 

IF (MOD (I— 1 , JBLKSZ)  . EQ .  0 ) RESMN ( I ) =RESMN ( I )* JBLKSZ 
DO  200  J=1 , N 
200  R ( J)  =  XMEAN(I) 

DO  250  J=1 , NSF 
S=XV( J) 

RT2=0 . 0 


=  1,N1 


NOW  GET  XMITTER  RESIDUAL 
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***  Subroutine  for  ADPCM  /  Kalman  filter  /  Zeroing  *** 
DO  230  K=1,N 

230  RT2=RT2  +  CM(K,I)*R(K) 

E2  =  S  -  RT2  -  RESMN(I) 

EQ2=QNTZ (E2 , NLEVEL , QL, QTH, VAREST ( I ) ) 

IF ( IB EAT  .NE.  1)EQ2=0. 

IF  ( IB  EAT  .EQ.  LIMCNT)  IBEAT*0 
IB  EAT* IB  EAT+1 
XTR  ( I ,  J)  =  EQ2 


SHIFT  R-REGISTER 


DO  240  K=2,N 
240  R (N+2-K)  =  R(N+1-K) 

R ( 1)  =  RT2  +EQ2  +RESMN { I ) 

250  CONTINUE 
160  CONTINUE 

VMN=VSME/256. 

WAR=(VSMSE-(  (VSME*VSME)/256.)  )/255. 

RMN=RSME/256 . 

RVAR= (RSMSE- ( (RSME*RSME) /256 . ) ) /255 . 

WRITE  (6,920)  VMN,  WAR,  RMN,  RVAR 

920  FORMAT  ( 1  VMN* 1  ,  E12 . 5  ,  ‘  WAR*  *  ,  E12 .5 , 1  RMN* ' , E12 . 5 , ’ 
WRITE (6 , 921) VLO, VHIGH 

921  FORMAT ( '  VLO* ' , E15.5 , '  VHIGH* ' , E15. 5 ) 


RECEIVER  SIMULATION  SECTION 


SET  R  -  REGISTER  TO  MEAN  VALUE 
FOR  THE  CURRENT  LINE. 


DO  210  I  *  1 , NF 
DO  170  J  =  1 ,  N 
170  R ( J) =XMEAN ( I ) 

DO  400  J  =  1 , NSF 
S  =  XIN ( I , J) 

RV2  =  0.0 
DO  180  K  =  1 ,  N 
180  RV2  *  RV2  +  CM(K,  I)  *R(K) 

RN  =  XTR ( I , J)  +  RV2  +  RESMN(I) 
XIN  ( I ,  J)  =  RN 


SHIFT  R  -  REGISTER. 


* 


x. 


RVAR*' ,E12.5) 

•*"«  ' 
i  ,  * 

*■  .  ’ 

•-'■I 


rj 
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***  Subroutine  for  ADPCM  /  Kalman  filter  /  Zeroing 


*  *  * 


DO  190  K  =  2 ,  N 
190  R (N+2-K)  =  R (N+l-K) 
R  ( 1)  =  RN 
400  CONTINUE 
210  CONTINUE 
RETURN 
END 


ooo  n n n o n n n o n n n n o o n n o n n n n o o 


***  Subroutine  for  ADPCM  /  Kalman  filter  /  Back-Quantization  *** 

SUBROUTINE  COMPMC (XIN, NF, NSF, N, CM, QTH, QL, NLEVEL) 

REAL  A  (10)  ,  VI  (10)  ,  XIN  (256, 256)  ,  R  ( 10 )  ,  XV  ( 256 )  ,  G  ( 10 )  ,  V  ( 10 , 10 ) 
REAL  XTR (256 , 256 ) 

REAL  CM ( 7 , 2  5 6 ) , QL ( 1 ) , QTH ( 1 ) ,E2V(256)  ,XMEAN(256) 

REAL  RESMN ( 256 ) 

REAL  QML (1024) , QMTH ( 1024) 

REAL  B ( 10 , 10 ) , PAR (10) 


EXPLANATION  OF  CALL  VARIABLES: 

XIN  -  THE  INPUT  DATA  MATRIX  AND  RESIDUAL  OUTPUT  MATRIX 
NF  -  NUMBER  OF  FRAMES  (NUMBER  OF  LINES  IN  THE  INPUT  MATRIX) 

NSF  -  NUMBER  OF  SAMPLES  PER  FRAME  (SAMPLES  PER  LINE) 

N  -  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED  ; 

CM  -  MATRIX  THAT  CONTAINS  THE  PREDICTOR  COEFFICIENT  VECTORS 

DEFINITION  OF  VARIABLE  TERMS: 

V  -  THE  ERROR  COVARIANCE  MATRIX 

A  -  THE  PREDICTOR  COEFFICIENT  VECTOR 

VARI  -  INITIAL  VALUE  FOR  THE  ERROR  COVARIANCE  MATRIX  ; 

W  -  VARIANCE  OFFSET 

XV  -  THE  INPUT  LINE  TEMPORARY  VECTOR 

G  -  THE  GAIN  VECTOR 

R  -  THE  PAST  VALUE  VECTOR 

E  -  THE  ERROR  OR  RESIDUAL  TERM 


SET  UP  THE  CONSTANTS  AND  INITIAL  VALUES 

DATA  V/100*0.0/,A/1.0,-.5,-.2, .3, .4, -.5/ 

CALL  MAXQTZ ( QMTH , QML , ML EV EL , 6 ) 

VLO=9999999. 

VHIGH=-1. 

XNSF  =  NSF 
COMMON  /A/  JBLKSZ 
XAVBLK= JBLKS Z  *1 2  8 . 

NSFMl=NSF-l 

W  =  1.0 

IN  THE  DO  160  LOOP,  I  IS  THE  LINE  NUMBER  (1  - 

DO  160  1=1, NF 
VARI=100. 

DO  10  J=1 , N 
DO  10  K=1 , N 
V ( J, K)  =  0.0 

IF(J.EQ.K)  V(J,K)  =  VARI 
CONTINUE 


non  nnnn  nn 


***  Subroutine  for  ADPCM  /  Kalman  filter  /  Back-Quantization  *** 

SET  UP  THE  INPUT  VECTOR  AND  THE  PAST  VALUE  VECTOR 

SUMX=0 . 

DO  20  J=1 , NSF 
SUMX=SUMX+XIN(I,  J) 

20  XV<J)  =  XIN ( I / J) 

XMEAN ( I ) =SUMX/XNSF 
WRITE (6,2000)1, XMEAN ( I ) 

2000  FORMAT ( '  LINE:  ' , 15 , 2X,  El 5 . 5 ) 

IF (MOD ( 1-1 , JBLKSZ )  .NE.  0)GO  TO  373 
XMEAN (I) =2750. 

GO  TO  374 

373  XMEAN ( I ) =0 . 

374  DO  30  J=1 , N 
30  R( J)  =  XMEAN ( I ) 

IDENTIFICATION  LOOP  (IDENTIFY  THE  PREDICTOR  COEFFICIENTS) 


DO  100  J=1 , NSF 
S  =  XV  (J) 

S2  =  0.0 
RT  =  0.0 
DO  40  K=1 , N 
VI (K)  =  0. 

DO  40  L=1,N 

40  VI (K)  =  VI (K)  +  V(K,L) *R(L) 

DO  50  K=1 , N 

50  S2  =  S2  +  R(K)  *V1  (K) 

DO  60  K=1 , N 

G  (K)  =  VI  ( K )  /  ( VAR  I  +  S2) 

60  RT  =  RT  +  A(K)*R(K) 

DO  70  K=1 ,  N 
DO  70  L=1,K 

V(K, L)  =  V(K, L)  -  G (K) *Vl { L) 

70  V(L,K)  =  V(K,  L) 

E  =  S  -  RT 
DO  80  K=1 ,  N 

80  A (K)  =  A(K)  +  G(K)*E 

SHIFT  THE  PAST  VALUE  VECTOR 
DO  90  L=2  ,  N 

90  R (N+2-L)  =  R (N+l-L) 

R ( 1)  =  RT+E 
100  CONTINUE 

CALL  PARSB L ( A , N , PAR , B ) 

DO  335  KK=1 , N 

CM (KK, I)=UNIFRM(-1.75,1.75,256,A(KK) ) 
335  CONTINUE 


lb  8 


***  Subroutine  for  ADPCM  /  Kalman  filter  /  Back-Quantization 


C 

C  RESET  R-REGISTER  TO  INITIAL  SIGNAL  VALUE 

C 

DO  120  J=1 , N 
120  R(J)  =  XMEAN(I) 

DO  150  J=1 , NSF 
S  =  XV  (J) 

RT2  =0.0 
C 

C  DETERMINE  RESIDUAL  SIGNAL 

C 

C  RT2  IS  AT  THE  TRANSMITTER  STAGE  WHERE: 

C  RT2  =  RT2  +  AO  (K)  *RV  (K) 

C 

DO  130  K=1 , N 

130  RT2  =  RT2  +  CM(K,I)*R(K) 

E2V(J) =S-RT2 
C 

C  SHIFT  R-REGISTER 

C 

DO  140  K=2,N 
140  R (N+2-K)  =  R(N+1-K) 

R ( 1 )  =  RT2+E2V(J) 

150  CONTINUE 
SUME=0 . 

SUMSQE=0 . 

DO  155  M=1 , NSF 
SUME=SUME+E2V(M) 

155  SUMSQE=SUMSQE+(E2V(M) *E2V(M) ) 

RESMN ( I ) =SUME/XNSF 

IF (MOD ( 1-1 , JBLKSZ )  . EQ .  0) RESMN ( I ) =RESMN ( I ) /JBLKSZ 

WRITE (6 »  91 9) I , RESMN ( I } , (A(K1) ,K1=1,N) , (CM(K2, I) ,K2=1,N) 

919  FORMAT (  1  '  ,  14 , 4  ( IX,  E10 . 3 )  ,/ ,  3  ( IX ,  E10 . 3 )  ) 

RSME=RSME+RESMN ( I ) 

RSMSE=RSMSE+RESMN  (I)  *RESMN(I) 

RESMN ( I ) =QNTZ (RESMN (I) , MLEVEL , QML ,  QMTH ,150.) 

IF (MOD (1-1, JBLKSZ)  . EQ .  0 ) RESMN ( I ) =RESMN ( I ) *JBLKSZ 
DO  200  J=1 , N 
200  R(J)  =  XMEAN(I) 

DELTA=100 . 

ALPHA1  = .  7 
AL  PH  A2  =  1 . 

RESMSQ=0 . 

DO  250  J=1 , NSF 
S=XV(J) 

RT2=0 .0 
C 

C  NOW  GET  XMITTER  RESIDUAL 
C 

f  DO  230  K=1 , N  1 


nnnooonono  ooooo 


***  Subroutine  for  ADPCM  /  Kalman  filter  /  Back-Quantization 


*** 


230  RT2  =  RT2  +  CM(K,I)*R(K) 

E2  =  S  -  RT2  -  RESMN(I) 

EQ2=QNTZ (E2,NLEVEL,QL,QTH, DELTA) 

XTR ( I ,  J) =EQ2 
E2V(  J)  =EQ2 
DELTA® 0 . 

DO  701  JI=1,J 

701  DELTA®DELTA+ (ALPHAl** (JI-1) ) *ABS (E2V ( J-JI+1) ) 
DELTA® ( ( 1-ALPHA1 ) /ALPHA2 ) * DELTA 
DELTA=DELTA*DELTA 
IF (DELTA  .LT.  VLO) VLO® DELTA 
IF (DELTA  .GT.  VHIGH)  VHIGH=DELTA 


SHIFT  R- REG IS TER 


DO  240  K®2 , N 
240  R (N+2-K)  =  R (N+l— K) 

R (1)  ®  RT2+EQ2+RESMN (I) 

250  CONTINUE 
160  CONTINUE 

VMN=VSME/256 . 

WAR=  ( VSMSE-  (  ( VSME*VSME)  /256  . )  )  /255 . 

RMN=RS ME/256. 

RVAR= (RSMSE- ( (RSME*RSME) /256 . ) ) /255 . 

WRITE  (6,920)  VMN,  WAR,  RMN,  RVAR 

920  FORMAT  ('  VMN®  '  ,  E12 . 5  ,  '  WAR®  '  ,  E12 .5  ,  '  RMN®  '  ,  El  2 .5  ,  '  RVAR=',E12 
WRITE (6,921) VLO, VHIGH 

921  FORMAT ('  VLO® ' , El 5 .5 , '  VHIGH® El 5. 5) 


RECEIVER  SIMULATION  SECTION 


SET  R  -  REGISTER  TO  MEAN  VALUE 
FOR  THE  CURRENT  LINE. 


DO  210  I  =  1 , NF 
DO  170  J  =  1 , N 
170  R ( J) =XMEAN ( I ) 

DO  400  J  =  1 , NSF 
S  =  XIN ( I , J) 

RV2  =  0.0 
DO  180  K  =  1 ,  N 
180  RV2  =  RV2  +  CM (K,  I)  *R (K ) 

RN  =  XTR ( I , J)  +  RV2  +  RESMN(I) 
XIN  ( I ,  J)  =  RN 


170 
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VAH PC  MAIN  DRIVER  PROGRAM 

EXPLANATIONS  OF  VARIABLES: 

XIN  -  MAIN  I/O  IMAGE  MATRIX 

COEFF  -  COEFFICIENT  STORAGE  MATRIX 

JBLKSZ  -  TRANSFORMATION  BLOCKS IZE 

ISZ  -  IMAGE  ROW  SIZE 

JSZ  -  IMAGE  COLUMN  SIZE 

NUMSMP  -  NUMBER  SAMPLES/LINE  (SAME  AS  JSZ) 

NUMFR  -  NUMBER  SAMPLE  LINES/FRAME  (SAME  AS  ISZ) 
LVECTR  -  NUMBER  OF  VECTORS  IN  CODE  BOOK 
KDIM  -  VECTOR  DIMENSION 

NCOEFF  -  NUMBER  OF  PREDICTOR  COEFFICIENTS 

SUBROUTINES: 

DSKIO  -  DISK  INPUT  /  OUTPUT 
COMPV  -  AD PCM  DATA  COMPRESSION 

USING  PARALLEL  PROCESSING 
AND  VECTOR  QUANTIZATION 


NUMFR=256 
NUMSMP=256 
NCOEFF =3 
LVECTR=I6 

LVECSQ=LVECTR*LVECTR 

LVECS1=LVECSQ-1 

KDIM=8 

ISZ=256 

JSZ=256 

JBLKSZ =16 

COMMON  /A/  JBLKSZ 

DO  10  1=1, ISZ 

CALL  DSKIO(A,JSZ,I,l,l,RL) 

DO  10  J=1,JSZ 

10  XIN ( I ,  J) =A( J) 

DO  11  1=1,12 

11  WRITE (6,1001) (XIN ( I , J) , J=1 ,16) 

1001  FORMAT ('  ' ,16 (F6.0,1X) ) 

CALL  COMPV (XIN, NUMFR, NUMSMP, NCOEFF, COEFF, T  VECTR, KDIM) 
WRITE (6,20) 

20  FORMAT ( '  ******  MADE  IT  BACK  THROUGH  SUBS  ******' 

DO  40  1=1, ISZ 


*  *  * 


Main  Driver  Program  for  VASPC 


★  ★  ★ 


DO  35  J=1,JSZ 
35  A(J)=XIN(I, J) 

CALL  DSKXO(A,JSZ,I,0,2/RL) 
40  CONTINUE 
WRITE (6, 2) 

2  FORMAT ('0***  JOB  FINISHED 
STOP 
END 


***  i  j 
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*** 


Max-Lloyd  Optimal  Quantizer  Design  Subroutine 


*** 


SUBROUTINE  MAXQTZ (QTH, QL, N,  NBIT) 


OPTIMAL  QUANTIZER  SUBROUTINE 

THIS  PROGRAM  CALCULATES  THE  OPTIMUM  QUANTIZATION 
LEVELS  FOR  A  SPECIFIC  NUMBER  OF  LEVELS  AS  OUT¬ 
LINED  BY  JOEL  MAX. 

THE  PROBABILITY  DENSITY  FUNCTION  IS  ASSUMED  SYM¬ 
METRIC  AND  IS  DEFINED  AS  A  FUNCTION,  HERE  AS 
FUNC. NORMAL,  AND  IS  ACCESSED  THROUGH  "FUNC". 


REAL*  8  X  (1024)  ,Y(1024)  , AREA, XLAST, THELST, XB ( 1024)  ,YB(1024) 

REAL *4  QTH(l) , QL ( 1 ) , ERROR 

REAL *8  EXTRA, ADD 

N=2**NBIT 

WRITE (6, 3) N 

3  FORMAT ('0NUM  LEVELS 15 ,/) 


THIS  SECTION  ZEROES  OUT  ALL  ARRAYS  AND  SETS  UP 
ALL  INITIAL  CRITERION  ESTIMATES. 


DO  5  K=1 ,1024 
X (K) =0  . 

XB  ( K )  =0  . 

YB  (K)  =0  . 

5  Y  (K)  =0  . 

XL AS T=999. 

ADD= .  1 

IF (N  .GT.  7 ) ADD= . 0 1 
IF (N  .GT.  20) ADD=,0005 
IF (N  .GT.  128) ADD=5.D-8 
GUESS=0 .0001 
TYPE=0 . 


THIS  SECTION  FIGURES  OUT  IF  N  IS  ODD  OR  EVEN. 

TYPE=1 .  — >  ODD 

TYPE=0 .  -->  EVEN 


CHECK=N/2 . 

ISTOP=N/2 

IF ( (CHECK-FLOAT (INT (CHECK) ) )  .NE.  0.)TYPE=1. 
IF (TYPE  .NE.  1 . ) GO  TO  10 
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***  Max-Lloyd  Optimal  Quantizer  Design  Subroutine  *** 

C - 

c 

C  THIS  SECTION  SETS  UP  INITIAL  ESTIMATES  OF  OPTIMUM 
C  SIGNAL  LEVELS  (X) ,  AND  QUANTIZATION  LEVELS  (Y) . 

C 

c 

Y(l) =0 . 

X (1) =GUESS 
THELST=X (1) 

Y (2) =2*X (1) -Y (1) 

GO  TO  20 
10  X(1)=0. 

Y(l) =GUESS 
THELST»Y(1) 

C - 

C 

C  WE  NOW  CALCULATE  ALL  SUCCEEDING  SIGNAL  AND 
C  QUANTIZATION  LEVELS  AS  OUTLINED  BY  MAX. 

C 

C 

20  DO  50  J=1 , ISTOP 
I=J+1 
JJ=J 

IF (TYPE  .EQ.  l.)JJ=J+l 
18  CALL  CENTRD(X(J) ,X(I) ,Y(JJ) ,AREA,0.) 

50  Y  ( JJ+1)  =2*X  (T)  -Y  ( JJ) 

C - 

C 

C  WE  NOW  SEE  IF  THE  LAST  Y-LEVEL  (QUANTIZATION)  IS 
C  THE  CENTROID  BETWEEN  THE  PREVIOUS  X-LEVEL 

C  (SIGNAL)  AND  INFINITY. 

C 

C 

JTRAN=JJ+1 

IEND=ISTOP 

IF (TYPE  .EQ.  1) IEND=ISTOP+l 

IF (N  .LE.  128) CALL  CENTRD (X ( ISTOP) , 15 ., Y ( IEND) , AREA, 1 . ) 

IF (N  .GT.  128) CALL  CENTRD(X ( ISTOP) ,100 ., Y (IEND) , AREA, 1 . ) 

IF (N  .GT.  128)  GO  TO  55 

IF ( (ADD  .LE.  l.D-6)  .OR.  (DABS (AREA)  .LE.  l.D-6))GO  TO  80 
GO  TO  56 

55  IF ( (ADD  .LE.  l.D-10)  .OR.  (DABS (AREA)  .LE.  l.D-6) ) GO  TO  80 

C - 

C 

C  IS  THE  /AREA/  WITH  THESE  VALUES  LESS  THAN  THE 
C  LAST  ITERATION? 

C 

C  IF  SO,  CHECK  TO  SEE  HOW  CLOSE  AREA  IS  TO 
C  ZERO,  AND  IF  NECESSARY,  MODIFY  INITIAL 

C  ESTIMATES  AND  ITERATE  AGAIN. 


nonnnnn  oonnn 


***  Max-Lloyd  Optimal  Quantizer  Design  Subroutine 


IF  NOT,  RESET  TO  LAST  ITERATION,  AND  REDUCE 
THE  MODIFICATION  VARIABLE  (ADD) . 


56  IF (DABS (AREA)  .LE.  DABS (XLAST) ) GO  TO  60 
IF (TYPE  .NE.  1.) GO  TO  57 

X (1) =THELST 

IF (XLAST  .EQ.  999.) GO  TO  58 
ADD2 ADD/ 2 . 

GO  TO  58 

57  Y (1) =THELST 

IF (XLAST  .EQ.  999.) GO  TO  58 
ADD=ADD/2. 

58  DO  59  J=1 , JTRAN 
X ( J) =XB ( J) 

59  Y ( J)  =YB  (J) 

60  XLAST2 AREA 

IF (TYPE  .NE.  1.) GO  TO  65 
THELST=X(1) 

GO  TO  67 
65  THELST2Y(1) 

67  IF (AREA  .LT.  0 . ) GO  TO  68 
EXTRA2 ADD 

GO  TO  69 

68  EXTRA2 -ADD 

69  DO  70  J=l, JTRAN 
XB ( J) =X(J) 

70  YB ( J) =Y ( J) 

IF (TYPE  .NE.  1.) GO  TO  75 
X ( 1 ) =X ( 1 ) +EXTRA 
Y (2) =2*X (1) 

GO  TO  20 

75  Y ( 1 )  =  Y ( 1 )  +i_ KTRA 
GO  TO  20 


WE  NOW  WRITE  OUT  THE  CALCULATED  OUTPUT  LEVELS 
—  SIGNAL  LEVELS  AND  CORRESPONDING  QUANTIZER 
LEVELS. 


80  NQL=N 
NQT=N-1 

IF  (TYPE  .NE.  1 . ) GO  TO  300 
NX=INT (N/2 . ) 

NY2 INT ( N/2 . )  +1 
FLG=0  . 

IPT=NX 

DO  205  1=1, NQT 


***  Max-Lloyd  Optimal  Quantizer  Design  Subroutine 


XXX=-X ( IPT) 

IF  (FLG  .EQ.  1 . ) GO  TO  201 
QTH ( I ) =-*X (IPT) 

IF  (FLG  .EQ.  0.)GO  TO  202 
1=1+1 

201  QTH  { I )  =X  ( I PT ) 

IPT=IPT+1 

GO  TO  205 

202  IF  (IPT  .EQ.  1 ) FLG=1 . 

IF  (IPT  .NE.  1) IPT=IPT-1 

205  CONTINUE 
FLG=0  . 

IPT=NY 

DO  210  1=1, NQL 
IF  (IPT  .EQ.  1 ) FLG=1 . 

IF  (FLG  .EQ.  1 . ) GO  TO  206 

QL ( I ) =-Y (IPT) 

YYY=-Y ( IPT) 

IPT= IPT-1 
GO  TO  210 

206  QL ( I ) =Y ( IPT) 

IPT=IPT+1 

210  CONTINUE 
GO  TO  999 

300  NX=INT(N/2 . ) 

NY=INT(N/2. ) 

FLG=0  . 

IPT=NX 

DO  305  1=1, NQT 
IF  (IPT  .EQ.  1) FLG=1 . 

IF  (FLG  .EQ.  1 . ) GO  TO  301 
QTH (I) =-X (IPT) 

IPT=IPT-1 
GO  TO  305 

301  QTH (I) =X (IPT) 

IPT= IPT+1 

305  CONTINUE 
FLG=0 . 

IPT=NY 

DO  310  1=1, NQL 

IF  (FLG  .EQ.  1 . ) GO  TO  306 

QL ( I ) =-Y ( IPT) 

IF  (FLG  .EQ.  0.)GO  TO  307 
1=1+1 

306  QL(I)=Y(IPT) 

IPT= IPT+1 
GO  TO  310 

307  IF  (IPT  .EQ.  1 ) FLG=1 . 

IF  (IPT  .NE.  1) IPT=IPT-1 
310  CONTINUE 
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FUNCTION  FUNC  (X ,  Q) 


NORMAL  GAUSSIAN  DENSITY  FUNCTION 


REAL* 8  FUNC, X, Q, P, ARG 
ARG--X*X/2 . 

IF (ARG  .LT.  -150.) GO  TO  1 
P=. 3 989422 804*DEXP (ARG) 

GO  TO  2 

1  P=0. 

2  FUNC= (X-Q) *P 
RETURN 

END 


FUNCTION  SQERR(X,Q) 


MEAN  SQUARED  ERROR  FUNCTION  FOR 
THE  NORMAL  GAUSSIAN  DENSITY. 


REAL *4  SQERR, X, Q, P 

P=.3989422804*EXP(-X*X/2.) 

SQERR= (X-Q) * (X-Q) *P 

RETURN 

END 
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***  Function  for  Adaptive  Quantization  *** 
FUNCTION  QNTZ (X, L, QL, QTH , VAR) 


ADAPTIVE  QUANTIZER  FUNCTION  WHICH 
QUANTIZES  INPUT  SAMPLES  GIVEN  THE 
DESIRED  LEVELS  AND  THRESHOLDS  ALONG 
WITH  THE  VARIANCE  ESTIMATE. 


REAL  QL (1) , QTH ( 1 ) 

LM1=L-1 
STD=SQRT ( VAR) 

DO  10  K=1 , LM1 
THRESH=QTH (K) *STD 
IF (X  .GE.  THRESH) GO  TO  10 
QNTZ=QL(K) *STD 
RETURN 
10  CONTINUE 

QNTZ=QL(L) *STD 

RETURN 

END 


**  Function  for  Uniform  Quantization 


#  *  * 


FUNCTION  UNIFRM (SMALL, BIG, N,X) 

IF (X  .LE.  SHALL) GO  TO  20 
IF (X  .GE.  BIG) GO  TO  40 
NM1=N-1 

ADD= { BIG-SMALL) /NM1 
TH  RADD= ADD/ 2 . 

THRSH2=SMALL 
DO  10  1=1 , NM1 
THRSH1=THRSH2 
DECIDE=THRSH1+THRADD 
THRSH2=THRSH1+ADD 

IF( (X.GT.THRSH1)  .AND.  (X. LE. DECIDE) ) GO  TO  11 
IF( (X.GT. DECIDE)  .AND.  (X. LE. THRSH2) ) GO  TO  12 

10  CONTINUE 
WRITE (6,1) 

1  FORMAT ( '  *.*.*.*>>  UNIFORM  LOGIC  ERROR  <<*.*.* 

11  UNIFRM=THRSH1 
GO  TO  99 

12  UNIFRM=THRSH2 
GO  TO  99 

20  UNIFRM=SMALL 
GO  TO  99 
40  UNIFRM=BIG 
99  RETURN 
END 
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***  Subroutine  for  Centroid  Calculations  *** 


SUBROUTINE  CENTRD (XI , XIP1 , Y, AREA, FLAG) 


CENTRD  IS  A  SUBROUTINE  DESIGNED  TO: 

1)  CALCULATE  THE  CENTROID  PIVOT  ELEMENT 
BETWEEN  TWO  ENDPOINTS  OF  A  KNOWN 
PROBABILITY  DENSITY  FUNCTION. 

2)  CALCULATE  THE  AREA  WITH  Y  AS  A  SPECIFIC 
CENTROID  PIVOT  ELEMENT. 


XI,  XIP1  ARE  THE  INTEGRAL  ENDPOINTS. 
Y  IS  THE  PIVOTAL  ELEMENT. 

-  XIP1 
I 

AREA* I  (X-Y) *P (X)  DX 

I 

XI 

FUNCTION  1)  IS  SELECTED  IF  FLAGOl. 
FUNCTION  2)  IS  SELECTED  IF  FLAG=1 . 


REAL* 8  XI , XIP1 , Y, AREA, YL, XLAST, ADTNL 
N=0 

IF (FLAG  .NE.  1 . ) GO  TO  5 
CALL  INTGRL (XI , XIP1 , Y, AREA, 0 . ) 

GO  TO  30 
5  XLAST=999. 

ADTNL*. 001 
CRTRN* .001 
XIP1=2. *Y-XI 
XL=XIP1 

10  CALL  INTGRL (XI, XIP1,Y, AREA, 0.) 

N=N+1 

IF (N  .EQ.  500 ) GO  TO  30 

IF (CRTRN  .LT.  DABS (AREA) ) GO  TO  1 5 

CRTRN* CRTRN/10. 

IF (CRTRN  .LE.  l.E-7)GO  TO  30 
15  IF (DABS (AREA)  .LE.  DABS (XLAST) ) GO  TO  20 
XIP1=XL 
ADTNL* ADTNL/ 2. 

CALL  INTGRL (XI, XIP1,Y, AREA, 0.) 

20  XLAST=AREA 
XL=XIPl 

IF (AREA  .LT.  0 . ) GO  TO  25 
X I PI =X I PI- ADTNL 
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***  Subroutine  for  Vector  Quantizer  Design 


it  it  it 


SUBROUTINE  VQDSN(Nf  KDIM, X, ISEQSZ ,A) 


VECTOR  QUANTIZER  DESIGN  SUBROUTINE 


A  =  QUANTIZER  CODEBOOK 
X  =  TRAINING  SET 

DSTRN  =  DISTORTION  BETWEEN  ITH  VECTOR  OF  X,  AND 
JTH  CODEVECTOR  OF  A. 

EPS  =  DISTORTION  THRESHOLD 

KDIM  *  VECTOR  SPACE  DIMENSION  (BLOCK  LENGTH) 

DM  =  OVERALL  AVERAGE  DISTORTION  AT  MTH  STAGE 
P  =  MINIMUM  DISTORTION  PARTITION 
N  *  NUMBER  OF  QUANTIZATION  LEVELS 
NUM  =  NUMBER  OF  TRAINING  VECTORS 
SMALL  =  LIST  OF  SMALLEST  DISTORTIONS 
IPT  *  IDENTIFIES  PARTITION  SET  FOR  EACH 
TRAINING  SET. 

NUMP  =  NUMBER  OF  TRAINING  VECTORS  IN 
EACH  PARTITION. 


DIMENSIONALITY: 

A ( N, KDIM) ,  X (KDIM, NUM) , 

DSTRN (NUM, N) ,  SMALL (NUM) ,  IPT(NUM), 
NUMP (N) 


REAL* 4  A(64,16)  ,X(KDIM,1)  ,SEQNCE(1) 

REAL* 4  SMALL (20 4 8) , ATEMP ( 32 , 1 6 ) ,EPSLON(I6) 
INTEGER  IPT(2048) ,NUMP(16) 

DATA  EPS LON/ 16 *.5/ 

DO  6  1*1,  N 
DO  6  K*1 , KDIM 
6  A { I , K) =0 . 

NUM= ISEQSZ /KDIM 
M=0 

EPS=0 .0001 
DLAST=1.E50 
DO  1  1=1, NUM 
DO  1  K=1 , KDIM 

1  A(1,K) =A(1,K) +X(K, I) /NUM 
NCOUNT=l 


SPLIT  CODEBOOK 
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***  Subcoutine  for  Vector  Quantizer  Design  *** 


5  DO  2  1=1 , NCOUNT 
DO  2  K=1 , KDIM 

2  ATEMP ( I ,  K) =A ( I , K } 

ICNT=0 

DO  4  1=1, NCOUNT 
ICNT= ICNT+1 
DO  3  K=1 , KDIM 

3  A(ICNT, K) =ATEMP {I , K) +EPSLON (K) 

ICNT= ICNT+1 

DO  4  K=1 , KDIM 

4  A ( ICNT, K) =ATEMP (I , K) -EPSLON (K) 
NCOUNT=2*NCOUNT 

10  IF  (M  .EQ.  100 )  GO  TO-  999 
DO  200  1=1, NUM 
SMALL ( I ) =1 . E5  0 
IPT ( I ) =1 

DO  100  J=l, NCOUNT 
DSTRN=0 . 

DO  50  K=1 , KDIM 

50  DSTRN=DSTRN+ (X (K, I) -A ( J, K) ) * (X (K, I) -A( J, K) ) 
IF (DSTRN  .GT.  SMALL (I)) GO  TO  100 
SMALL ( I ) =DSTRN 
IPT ( I ) =J 
100  CONTINUE 
200  CONTINUE 

DO  300  J=l, NCOUNT 
NUMP(J) =0 
DO  300  1=1, NUM 
I F ( I PT ( I )  .NE.  J)  GO  TO  300 
NUMP(J) =NUMP ( J • +1 
300  CONTINUE 
DM=0  . 

DO  400  1=1 , NUM 
400  DM=DM+SMALL { I ) 

DM=DM/NUM 
WRITE (6, 772) DM 
772  FORMAT ('  DM:', El  5. 5) 

IF (DM  .EQ.  DLAST)  GO  TO  999 

IF (ABS ( (DLAST-DM) /DM)  .LE.  EPS)  GO  TO  999 

DLAST=DM 

DO  500  1=1, NCOUNT 

IF ( NUMP ( I )  .EQ.  0)  GO  TO  500 

NV£CS=NUMP ( I ) 

DO  420  K=1 , KDIM 
420  A ( I , K) =0 . 

IF (NVECS  .EQ.  0) GO  TO  430 
DO  430  J=1 , NUM 
I F ( I PT ( J )  .NE.  I ) GO  TO  430 
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Subroutine  for  Vector  Quantizer  Design 


★  *  it 


DO  425  K=1 , KDIM 

425  A(I,K)=A(I,K) +X(K, J)/NVECS 
430  CONTINUE 
500  CONTINUE 
M=M+1 
GO  TO  10 

999  IF (NCOUNT  .NE.  N) GO  TO  5 
WRITE  (6, 888)  M 
888  FORMAT ( '  M» ' ,15) 

RETURN 

END 


r-'  • 
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*  *  * 


Subroutine  for  Vector  Quantization  *** 
SUBROUTINE  VECQNT(X, KDIM, N, CDBK) 


VECTOR  QUANTIZER: 

EXHAUSTIVE  SEARCH  FOR  MIN  DISTORTION. 


REAL *4  X(l) ,CDBK(64,16) 

SMALL-1 .  E50 
IPT=1 

DO  50  J=1 , N 
DSTRN=0 . 

DO  25  K=1 , KDIM 

25  DSTRN=DSTRN+ (X (K) -CDBK ( J,  K) ) * (X (K) -CDBK ( J, K) ) 
IF (DSTRN  .GT.  SMALL) GO  TO  50 
SMALL® DSTRN 
IPT=J 

50  CONTINUE 

DO  75  K»1 , KDIM 
75  X (K) »CDBK (IPT, K) 

RETURN 

END 
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**  ProyL'usr.  for  Tape  /  Disk  Transferals  *** 


INTEGER*2  IP(256) 

LOG ICAL*1  PPICT (2,256) 

EQUIVALENCE ( IP (1) ,PPICT<1,1) ) 

DO  999  N-1,6 
DO  10  1=1,256 
READ{N) IP 
DO  7  KK=1 ,256 

IF  (IP(KK)  .GT.  255)  IP(KK)=255 
IF  (IP(KK)  .LT.  0)  IP (KK) =0 
7  CONTINUE 

WRITE (8,11) (PPICT (2, MM) ,MM=1,256) 
11  FORMAT ( 4 ( 6 4A1 ) ) 

10  CONTINUE 
ENDFILE  8 
999  CONTINUE 
REWIND  8 
STOP 
END 
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***  Program  for  SNR  Calculations  and  Histogram  Generation  *** 
INTEGER*2  R ( 256 ) , T ( 256 ) 

INTEGERM  RI{256)  ,TI(256)  ,  ORIG  ( 256 , 256 )  ,  TWOISZ  ,  TWOIP1 
DATA  R/256*0/ 

REAL  H(12,513) ,TMSE(12) ,SNR(12) ,RATIO(12) ,HWR(257) 

REAL  SIG ( 16 ) 

DATA  SIG/16*0 ./ 

C - 

C 

c 

c 

PROGRAM  TO  CALCULATE  THE  SNR  RATIO  OF  AN  IMAGE 
AND  GENERATE  A  HISTOGRAM  OF  THE  ERROR  BETWEEN 
TWO  IMAGES. 


ISZ=256 
READ (5,17) NUM 
17  FORMAT (12) 

SQISZ=FLOAT( ISZ ) *FLOAT(ISZ) 

TWOISZ=2*ISZ 
ISZP1=ISZ+1 
TWOIPl=TWOISZ+l 
INCRMT=512-ISZ 
SUM=0 . 

DO  10  1=1, ISZ 
READ ( 1 ) T 
DO  10  J=1,ISZ 
TI(J)=T(J) 

SUM=SUM+FLOAT(TI(J) ) 

10  ORIG ( I , J) =TI ( J) 

AVE=SUM/SQISZ 

SIG2=0. 

1=0 

IC=0 

DO  16  ID=1 , ISZ ,16 
IC=IC+1 
DO  15  IN=1 ,16 
1=1+1 

DO  15  J=l, ISZ 

15  SIG(IC)=SIG(IC) + (FLOAT (ORIG (I, J) ) -AVE) *( FLOAT (ORIG ( I , J) )-AVE) 

16  SIG2=SIG2+SIG(IC) 

DO  200  N=1 , NUM 
NUMOFF=N+7 

DO  20  J=l, TWOISZ 
20  H (N, J) =0 . 

ERR2=0. 

WRITE ( 6,2) N 

2  FORMAT ( '  READING  IMAGE  #',I3) 


vss 


***  Program  for  SNR  Calculations  and  Histogram  Generation 
1C=0 

DO  50  I0UT=1 , ISZ ,16 

IC=IC+1 

ERR1=0. 

DO  49  IN=1 ,16 
1  =  1+1 

READ (NUMOFF) R 
DO  49  J=l, ISZ 
RI ( J) =R( J) 

IDIFF=ORIG(I, J) -RI ( J) 

IF ( IDIFF  .LT.  (ORIG ( I , J) -255) ) IDIFF=ORIG (I , J) -255 
IF ( IDIFF  .GT.  ORIG(I, J) ) IDIFF=ORIG (I,J) 
INDEX=IDIFF+ISZPl 
H (N, INDEX) =H(N, INDEX) +1.0 

49  ERRl=ERRl+FLOAT( IDIFF) *FLOAT( IDIFF) 

ERR2=ERR2+ERRl 

BLKSQ=40 96 . **2. 

TMS1 =ERRl/BLKSQ 
SNR1=SIG ( IC) /ERR1 
RATI =10. *ALOGl 0 ( SNR1 ) 

WRITE (6. 400) IC, TMS1 , SNR1 , RATI 
400  FORMAT ('  BLK:  ' , 15 , 3 ( 2X, El 5 . 5 ) ) 

50  CONTINUE 

100  TMSE (N) =ERR2/ SQISZ 
SNR(N) =SIG2/ERR2 
200  RATIO ( N) =10. *ALOG10 ( SNR (N) ) 

DO  300  N=1 ,NUM 
DO  250  J=1 , TWO I Pi 

IF  (  ( J  .LT.  128)  .OR.  (J  .GT.  384)  )  GO  TO  250 

JPT=J-127 

HWR( JPT) =H (N, J) 

250  CONTINUE 

WRITE (24)  HWR 

WRITE ( 6 , 1 )  N, TMSE (N) ,SNR(N)  ,RATIO(N) 

1  FORMAT ( '  #',13,'  MSE:  ',£15.5,'  SNR:  ',E15.5,'  R 
300  CONTINUE 
STOP 
END 


RATIO:  ' , El  5 . 5 ) 
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***  Subroutine  to  Plot  Histograms  —  CALCOMP  *** 

REAL  HWR(257)  ,H(259)  ,X(259) 

IB=1 

NC=1 

CALL  PLOTS  ( IB, NC, 25) 

CALL  PLOT  (1.75, 1.0, -3) 

LG=257 

NUM=6 

DO  500  LOOP=l , NUM 
WRITE (6,11) 

11  FORMAT ( '  /*CHANGE  PAPER  THEN  ENTER:  1*) 

READ (5,12) ANS 

12  FORMAT (1A1) 

READ (24) HWR 
WRITE (25, 13) 

13  FORMAT ('  J  COPY  HO  PL1 • ) 

DO  10  1=1, LG 

10  H ( I ) =HWR ( I ) 

DO  50  1=1, LG 
50  X (I) =1-129 

X  (LG+1)  =-128.0 
X (LG+2) =32.0 
H  ( LG+1 )  =0 . 0 
H (LG+2) =2200.0 

CALL  AXIS  (0. ,0. , 'PIXEL  ERROR* , -11 , 8 . 0 , 0 . 0 , X (LG+1 ), X (LG+2) ) 
CALL  AXIS  (0. ,0. , 'FREQUENCY  ' , 10 , 5 . 0 , 90 . , H ( LG+1 ) , H (LG+2) ) 

DO  100  1=1, LG, 1 
XX=(X(I)-X(LG+1) ) /X (LG+2) 

YY=H ( I ) /H (LG+2 ) 

CALL  PLOT (XX, YY, 3 ) 

100  CALL  PLOT (XX, 0. ,2) 

CALL  SYMBOL(l. 5, 5. 0,0. 21, 'HISTOGRAM  OF  THE  ERROR ', 0 . 0 , 22) 
CALL  SYMBOL (2.55,4.75,0.21, ' FOR  THE  CITY' ,0 .0 ,12) 

CALL  SYMBOL (3. 075, 4. 5, 0.21, '256X256 ’ ,0.0,7) 

CALL  PLOT(12. 0,9. 0,999) 

500  CONTINUE 
STOP 
END 

SUBROUTINE  PSPCHR 
C  SUB  FOR  SPECIAL  CHARACTERS. 

RETURN 

END 
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